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INTRODUCTION 


Classical techniques for analysis and design of dynamic 
systems are largely restricted to cases in which on-ly one para- 
meter of the system is adjustable. As a consequence complex 

systems cannot be treated adequately with classical techniques. 

* 

Algebraic methods,, as developed in NASA CR-617 . are capable, of' 
treating systems in which two parameters . are adjustable, and 
thus permit analysis, and synthesis of systems which are too 
complex for treatment with classical methods. 

The treatment of algebraic methods presented in CR-617 
develops the fundamental theoretical basis for the . coefficient 
plane and parameter plane methods. It also applies these 
methods to basic problems such as stability analysis 1 , cascade 
compensation of systems,' and related topics. The applications 
indicated in CR-617 are rather elementary, i..e., the problems 

considered illustrated the procedures to be used but- were not 

* »• 

very complex? p'roblems,. ‘‘This report - is based on the findings of 

* ‘ . 1 r * * 4 - 1 

CR-617 1 , and extends the applications of the alcrebrain ’methods to 
problems of a mbre complex, nature,. 

When cascade .compensation is used in a feedback control 
system, more than one- filter section may be required to achieve 
desired performance. Frequency response methods involving trial 
and- error are often used, but parameter plane methods permit 
analysis and design without trial and error if it is permissible 
* 

Algebraic Methods for Dynamic Systems by G„ J.. Thaler, D. D. 
Siljak and R. C. Dorf, Nasa Contractor Report NASA CR-617, 

Nov., 1966. 




to use two identical filter sections. This problem ,is treated 
in Chapter I of this report. The applicable parameter plane 
equations are derived and a digital computer program based on 
these equations is presented. The program 1 is- used to study the 
effects of compensation on several systems. 

Chapters II and III are concerned with nonlinear systems. 
Conventional methods such as frequency domain analysis of sys- . 
terns with the Describing function have proven useful when the 
system contains only one nonlinearity (or several nonlinearities 
conveniently located so that they can be incorporated in one de- 
scribing function) . These techniques can define stability and 
estimate relative stability for fairly complex systems as long 
as the -conditions- of , nonlinearity are not 'too complex. Such 
cases are^easily treated using algebraic methods, the effect of 
the nonlinearity being represented as a movement -of the operating 

, - _ - ■ > , i . 

point on the parameter plane', which' in ’ turn represents a varia- 
tion of the characteristic, roots as a function of signal ampli- 
tude. The algebraic methods are capable of extending such analy- 
sis to systems containing two distinct nonlinear components, and- 
ean be used to predict the transient response of the system 
rather accurately. Techniques for such problems are developed 
in Chapter II. 

Chapter III is concerned with a much more difficult non- 
linear problem, that of asymmetrical nonlinear oscillations. 

These are oscillations consisting of a limit cycle superimposed 
on another signal. The problems studied on the parameter plane 



involve .steady-state operating conditions (rather than transient 
conditions) , and permit analysis of the existence of oscil- 
lations as well as their dependence on parameter values and in- 
put signal values. Extension to linearization with either 
signals is included, as well as some design considerations. 

It. is felt that the results presented here indicate the 
•capabilities of the algebraic methods in dealing with complex 
linear and nonlinear problems. It is also felt that the re- 
sults presented here will be .directly applicable to a number of 
practical problems, and will point out avenues of approach to 
still additional problems. 



I 


SOLUTION OF EQUATIONS WITH COEFFICIENTS 
THAT ARE QUADRATIC IN or and B 


1 . 1 INTRODUCTION 

It has been shown that the characteristic equation can be 
solved for or = and B = when the coefficients of 

the characteristic equation are of the forms: 


a> ®k “ V + c * s + ^ 

b) = b k « + cyS + rfi + ^ 

c) ®k = \2 a2 + \l a + V S + °*10 + c k2^ 2 + ^ 


( 1 ) 


d) *k = b kn 


a n + b a n 1 + 

k(n-l) 


'k(n-l) 


S "' 1 + + \ 




In addition practical solutions have been obtained for the first 

. ‘ ' ’ ■ . ' , i _ . , ■ 

two of these coefficient forms, i.e., computer programs have been 
written for them and successfully applied. The development to 
be presented here, is a particular solution for case 1-lc, parti- 
cularly in the sense that a computer program has been obtained 
which solves the equations of a third order system for which tbe 
coefficients are quadratic in or and 8, but which do not contain 
all of the o' and 8 ‘combinations indicated. At the same, time the 

solution is a general solution in the sense that the program ca 

• « til ♦ ^ 

be modified to solve the equations of an n order -system, aiR 

can also be modified to accept all of the oc and 8 forms^ indi- 


cated in 


** b k2^ + b kl ff + \. aB + c kl 8 + c k2 8 + \ 


k2 ~kl Tc kl k2 

The modifications to be made in- the program are discussed, but 
the necessary programming has not been done. 



1.2 THE PROBLEM: Cascade Compensation with two identical 

filter sections . 

In the design of feedback control systems it is common to 
use compensators which are filters placed in cascade with the 
main transmission path. Frequently two sections of filter are 
needed, and if identical sections are used with an isolation 
amplifier so that their transfer functions can be multiplied, 
then manipulation of the transfer function equation provides a 
characteristic equation in which the coefficients- are quadratic 
in z and p, the zero and pole of the compensators. For example 
let: 


G - 

G = 
c 


K 


0,3 2 

§ +Xs +Ys 


(: 


s + z )2 = 


s+p' 


2 2 

s +2zs+z 

2 - 2 

s +2ps+p 


(1-2) 

(1-3) 


1+G G = O' = 
c t 


, - K{s^+2zs+Zi 2 y)(f ■ - 

3 ,2 — : — 2-' 1 — T 

J ( s +Xs' +Ys )’ (s t2ps+p ) 


(1-4) 


from which the character! Stic equation is 

s 5 + (X+2p) s 4 + (p 2 +2Xp+Y) s 3 + , (Xp 2 +2Yp + K) s 2 + 

\ ’ ' - ‘ f ■> ’ . 

• + (Yp 2 + 2Kz)s + Kz 2 = 0 (1-5) 

. A A 

Letting p — o' and z = 8 it is noted that all of the forms speci- 
fied in the quadratic case definition of a^, do appear in at leasi 

some of the coefficients except that there is no &8 product term. 

The formulation just given does not conform to normal con- 
trol system practice, however, in that an important restriction 
on the design of the compensator is the usual requirement that 
steady state accuracy must be maintained by keeping the error 
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In equation 1-8 the coefficients are quadratic in & and 8 , 
but there is no term of the form h^o?, and the program as written 
does not make provision for such a term, though modification of 
the program to- include it is not difficult. The problem to be 
studied, then is that of a third order system compensated with 
two cascaded identical sections of filter, and with the added 
requirement that the error coefficient be maintained constant 
at a predetermined value. 



1.3 DERIVATION OF THE GENERAL THIRD ORDER SYSTEM 
RELATIONSHIPS 


The general third order system is defined by the transfer 
function 


G(s.) = 


K 


(s+A) (s+B) (s+C) 

which is a Type Eero system,, but which can be changed to Type 
1, 2, or 3 by setting one or more of the poles to zero. The 
compensator transfer function-, including the -gain multiplier 
which maintains the error coefficient is 


(1-9) 


2 2 2 
,p„ 2 ,s+z x 2 p (s +2zs+z ) 

G = (tr) ( — - — ) = x tt 

c z v s+p 2 , 2,. , 2. 

* z (.s +2ps+p ) 


( 1 - 10 ) 


From 1-9 and 1-10 the characteristic equation is 

[s 3 + (A+B+C)s 2 + (AB+BC+AC) s + ABc] (s 2 +2ps+p 2 ) + 


+ K ^ (s 2 +2zs+z 2 ) = 0 

V. 


d-11) 


This expands to- 


+ ( A+B+C+2p ) s 4 \ + T AB+BC +AC +2p ( A+B+C ) +p 2 Is 3 

* , i * 

* t 1 , * 1 * * ' ♦ *■ 2 

+ [abC+ 2 P (AB+BC +CA) + p 2 .(A+B+C) + K ^Js 2 

z 

+ [2pABC+p 2 (AB+BC+AC) + 2Kp 0|) ]s+p 2 (ABC+2K) = 0 


A 

Let p = $ 


E A 


= ' a 


A+B+C 

AB+BC+AC 


ABC 


— ^ r. = sum of roots (poles) 

A V rr 

- £ = sum of root products taken 2 at a time 

^ Or^ ^ sum of root products taken n at a time 

= ITr^ = products of the roots 
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Then equation 1-12 becomes: 

■s 5 + (£r ± + 2P)-s 4 + (£ H r ± + 2^r ± $+$ A ) s 3 + 

(H r.+2£ n r ± P + Jr ± $ 2 + K« 2 )s 2 + (1-13) 

(2n r ± P + Y, -H r ± P 2 + 2K<*P)s + (II r ± +2K) p 2 = 0 
Collecting like terms in or and P': 

o; 2 (Ks 2 ) + ap(.2Ks) + P 2 (s 3 + Yr.s 2 + Y nr.s + n r.+K) + 

L-j 1 i — t ^ J- 1 

+ p (2s 4 + 2^r i s 3 + 2^ .n r^s 2 + 2l3r^s) 

\ *« 

+ (s 5 + Yr.s 4 + Y Hr . s 3 + Tlr.s 2 ) = 0 (1-14) 

J—t 1 l—i ^ 1 1 


Using the basic parameter plane relationships : 
n 

^ = 0 

k=0 

n. 

£ (-l) k U k (C) = 0 

h=0 

and defining: 

B 21 " K * 2 °l(?) 

b 22 = k» 2 d 2 (4) 

D, = -2K®U (£), 

1 o 

D 2 = -^KWU^E)' 

E 1± = 2wV ( (£) -2^.-, w 3 U 2 (?) + 2^ n rVu^O 

“2nr i «U o C0 

e 12 = 2w 4 u 4 (a -2^10^3(0 + 2 £ n r ± «) 2 u 2 (s) 

-2n r ± wu l (^) 


I <- i) ' k ^ 


(1-15) 

(1-16) 

(1-17) 

(1-18) 

(1-19) 

1 - 20 ) 

1 - 21 ) 

. 1 - 22 ) 
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(1-23) 

(1-24) 

(1-25) 

(1-26) 

(1-27) 

(1-28) 

(1-29) 

(1-30) 

(1-31) 

(1-32) 


which are two non-linear algebraic equations completely general- 
ized in terms of the uncompensated system poles and root locus 
gain, £, and' the first kind of Chebyshev Functions. These 
must be solved simultaneously for the correct values of ot and 
P . To do this-, the method with the best chance of success 
appears- to be Sylvester's Method in which we form a set of four 
equations by taking the original Equations (1-31) and (1-32) and 


forming two more by a multiplication with ot giving: 



0 

0 


6 


(1-33) 

(1-34) 




a' 

B 21 

- " 2p i 

+ 

= 0 




(1-35) 


a' 

S 2 

+ a' 2 P 2 

+ <*q 2 

= 0 




(1-36) 

Now placing 

these 

equations in 

matrix form: 










r* ^ 





0 

B 21 

P 1 

Q l 


or 3 





0 - 

B 22 

F 2 

Q 2 : 


. 2 
■a 





B 21 

P 1 

Q 1 

0 


(X 

= 0 

(1-37) 



_ B 22 

1 p' 
2 

: Q 2 • • 

0 ’ 

■ 

1 

■ 



If the or ' s are not' zero : 

then: 

; t 

f 

t . 





0 

B 21 

P 1 

Q 1 







0 

B 22 

P 2 

Q 2 







B 21 

P 1 

Q 1 

0 

“ 

= 0 


C1-38-) 



B 22 

P 2 

q 2 ■ 

0 





Expanding this determinant 






2 2 
- B 21 Q 2 

+ 

B 21 B 22^1 Q 2 

+ p p 
12 

B 21 Q 2 - 

P l°2 

B 22 + Q 1 Q 2 B 21 B 22 






2 2 

- Q 1 B 22 

Q 1 P 2 B 21 + 

Q 1 P 1 P 2 B 22 = 0 

(1-39) 


Substituting equations (1-27) through (1-30) in equation (1-39) 


provides a fourth order equation in P : 

P ^" F 22 B 21 + 2F 21 F 22 B 21 B 22 + D 1 D 2 F 22 B 21 “ F 22 D 1 B 22 


2 2 
F 21 B 22 


F 21 D 2 B 21 + D 1 D 2 F 21 B 22 ) + 


g (-2E 12 F 22 B 21 + 2 E i;l F 22 B 21 B 22 + 2F 2l E l2 B 21 B 22 + 
D 1 D 2 E 12 B 21 _ D 1 E 12 B 22 “ 2 ’ E 11 F 21 B 22 


D 2 E 11 B 21 + D 1 D 2 E 11 B 22 ) + 

^ ^“ E 12 B 21 " 2F 22 G 2 B 21 + 2E 11 E 12 B 21 B 22 + 2F 21 G 2 B 21 B 22 + 
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2G 1 F 22 B 21 B 22 + D 1 D 2 G 2 B 21 _ D 1 G 2 B 22 “ E 11 B 22 “ 
2F 21 G 1 B 22 “ D 2 G 1 B 21 + D 1 D 2 G 1 B 22^ + 

P ( -2E 12 G 2 B 21 + 2E 11 G 2 B 21 B 22 + 2G 1 E 12 B 21 B 22 “ 2E 11 G 1 B 22 J + 


+ 2G I G 2 B 21 B 22 " G 1 B 22 ) = 0 


(1-40) 


from which the coefficients may he determined by a substitution 
of '(1-17) through (1-26) 'and 'the ''values 1 of the first kind of 
Chebyshev functions in i terms of- £ and w. -Since the solution of 

T « ' „ t * 

t . » - " 

a fourth order equation is at best difficult, it is at this 

. .... ' * * ' » i 1 

point a digital computer becomes a' necessity* 

The major problem is not the actual solution of the quartic 
itself, but rather the proper choice of one of the four solutions. 
There are two marked characteristics, however, which help in the 
selection. These are: 

a) Complex answers to the quartic have no physical signifi- 

J 

cance and may therefore be discarded as erroneous. 

b) The definition of a requires that o' and p be of the 

same sign so that p and z will be of identical sign. 

. . 2 
Using this information and that available from the Ross-Warren 

method as to compensator pole and zero location, if is found that 

the solution to the P quartic is the largest, positive, rea 

value. 


Now entering equation (1-27) with this value, and evaluating 


the other coefficients 


for in the third order case P n is always identically zero. 


(1-41) 
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Thus, with the programming of the appropriate equations, 
the digital computer could give all of the values and plot the 
constant zeta and constant omega loci on the Parameter Plane for 
any desired values. 

1.4 SOME APPLICATIONS OF THE. PROGRAM 

Several third order systems were investigated by the appli- 
cation of the generalized equations and the Parameter Plane 
curves. Figures 1-1 through 1-8 were plotted. Of these, the 

3 

K/s family appears the most interesting. Further investigation 
of three of the curves in this family. Figures 1-1, 1-2 and 1-3 
shows that there is a relationship between K, the root locus 
gain, a and P. 

These relations are: 

3 

a) Choose a point on the 1/s ct?-p plane. 

b) Zeta reads directly. 

c) Determine the actual omega at that point by multiplying 
the value read by the cube root of the uncompensated 
system gain. 

d) Read the value of ct directly from the point chosen. 

e) Read the value of P from the point chosen. 

f) Obtain the true value of P by multiplying this value 

by the cube root of the uncompensated system gain. 

By this method, the values of a and P may be determined for 
K 

all — ^ systems from one universal curve. 
s J 

1.5 BANDWIDTH CURVES ON THE Q'-p Plane 

In many instances, there is also a bandwidth criterion 
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imposed on the engineer as well as an optimal operating point for 
the plant under consideration. With this in mind, equations for 
the plotting of constant bandwidth curves on the Qf-P plane' are 
developed. For the purpose of this development a constant band- 
width curve will be defined as : 

A constant bandwidth curve for G(j^) = M is a cu-rve drawn 
upon the parameter plane which specifies the relation between 
the parameters necessary if the transfer function G(s) , which 
is a function of the parameters, is to have magnitude M at 
the real frequency 

Once these curves are obtained they may be superimposed on the 
parameter plane thus indicating what values of the parameters are 
necessary in order to meet the specifications. 


Taking the rational transfer function and defining it: 

m , m— 1 , _ _ , „ 

_ , . ps +P . S + ... Pi S- + p_ 

f ~ * r£ - ! P (s) _ f_m fjm-1 ^1 o 

q 'S 11 + q . s n ^ + . . . q n s + q ' 
H n ^n-1 "1 o 


-G-(-S-)' -= 


Q(s') 


(1-42) 


where the p 's and q 's are of the form: 


p = go , ^+h® + io'P+jP+kP 2 + l 

* u ^u u u u u u 


u = 0, 1, 2, ... ,m 

q - a a'i + b c +, d P + ;-e j $' 2 + f 

^v J v , 1 ■ V ‘ 1 V V , , , V . * V 


(1-43) 


Therefore 


t P, 


u 


g,<sj; 


u=0.’ 


1 

v=0 


V 


v — ■ 0,1 ,2 ,. . . , n 


iu 


(1-44) 


(1-45) 



Employing Equation 1-45 in the parameterized form the generalized 
compensated third order transfer function is: 


G ( s y = P.( s .).. 
G(S) Q (s)' 


(1-46) 


P(s) = Q' 2 Ks 2 + 2 c*PKs + P 2 K 


(lr47) 


and: 


Q<s) = P 2 [s 3 + l n r ± s 2 ^nrs+^B rj + 

x 2 3 

p |^2s^ + 2 n r_^s 3 + 2 ^ n r^s 2 + 2 ^ n r^sjn 

(s 5 + ^ n r^s^ + ^ II r^s 3 + ^ II r^s 2 J (1-46 


Making the definitions : 


= y (-l)^u^a ; etc. for B ,C ,D ,E ,F 
L> ' ' b v r r r r r 


(1-49.) 


even 


L i = I (-1)J5(V_1) “S a v ? etC * for V C i' D i' E i' F i (1-50) 




V ^ 1 

Gj: = Z (-l)^ U ^g u 7 etc. for 


(1-51) 


even 


ILL 

G. =7 (-l)'^ U-1) wjjg ? etc. for (1-52) 

l L, b^u i ilia. 


and substituting in Equation 1-46, 

(a 2 G + K ) + j (apl.) 

G(jV = — o — — 2^ 

O D + PE + F ) + j (PD ± + PE ± + F ± ) 


(1-53) 



Setting the magnitude of G(j<^) = M: 


(Q? 2 G + K ) 2 + (api.) 2 

. r r v x 


2 ^ 5 

(P.D. + PE + F ) ' + (P D. + PE. + F.) 

* "i a? r r' t. v x 1 x' 


Manipulatinq Equation (1-54) alqebraically 

0(or # P) - M 2 0 (a,p) = 0 (l-5£ 

where fy 2 2 2 222 

0(g-,P) = ff 'G + 2« K G + K7. + if (l-5€ 

P -L JL, JL JL 1 

0(<*,P) = P 4 D 2 + 2P 3 D E + 2p 2 D F + p 2 E 2 
r rr rrr 

2pE F + F 2 + P 4 D? + P 2 E? + F? 
r, r r x x x 

2P 3 E.D. + 2 J3 2 D . F . + B 2 E? + F? 

XX XX XX 

2P 3 E.D. + 2p 2 D.F. + 2PE.F. (1-5: 

xx xx xx 

Substituting Equations' (1-56) and (1-^57) in Equation (1t5'5) and 


defining : 


2 2 
P, = D + B. 
1 r x 


2D E. + 2E.D. 
r r xx 


2D F + E 2 + E? + 2D , F . 
r r r 1 x x 


2 2 
or IT 
x 


2E F + 2E.F. 
r r xx 


2 2 

F + Ff 
r x 

X 

4 2 2 2 

a G + 2 a K G + K 
r rrr 
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It follows that 

M 2 p n p 4 + M 2 Q r P 3 + (M 2 R- - R 2 )P 2 + M 2 V x P + 

(M 2 W 1 - W 2 ) = (1-65) 

Since the Parameter Plane for compensation purposes has already 
been determined it is now a matter of taking the computed a 
values and substituting them along with a constant value of 
omega and M into Equation (1-65) and then solving the P quartic. 
This has as its solution the largest, real and positive value of 
the four roots as before. 


1.5 EXTENSIONS TO HIGHER ORDER SYSTEMS 

Although the work presented to this point has been limited 
to third order systems and the program written for this specific 
case, investigation shows that generalized equations may be 
written which will allow the extension of:. 'the program to higher 
ordered systems. It can be shown for a given n order system 
with no zeros to be compensated with two identical sections of 
cascade compensation, that the characteristic equation of the 


system may be generally written as: 

n+2 , „ n+1 , 2n,,22, 0 , 2. 

s + 2ps + p s + (z s + 2pzs + p )K + 

j=n j=n 

k=l k=0 

2 b I (I n r i> s* + p 2 £ .(£ n r ± ) s* + 

k=n 2 3 

j=l 3 A 

k=2 

X <x ? r i> slc = 0 

k=n+l : 
j=l 
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where for n=4 the equation would be written: 

6,0 5 2 4 2 2 _ 2 . 

s + 2 ps + p s + (z s + 2 pzs + p- )K + 


2p (Y . II r . s 4 +Y ITr.s 3 + ') llr.s 2 + Y H r . s) 
L i 1 ^ 2 1 ^ 3 1 ^4 

p 2 <1 ? r ± s3 + 1 ? v 2 + 1 ? v + 1 ? r i> + 


i , 


n r^s 5 + y n r^s 4 + n r..s 3 + ^ ^ r^s 2 ) 


(1-67) 

It may be further shown that the parameters defined by Equations 


(1-17) through (1-26) may be written: 
B,, = k«j 2 u. ux. 


B 22 = K “ V«> 

D, = -2K«JU (£) 
1 o ^ 

D 2 = - 2 X ^(0 


11 - 68 ) 


(1-69) 


(1-70) 


(1-71) 


j'i. — 

= 2 (-l) n+ V +1 U n (0. + 2 I [.(-i^^iC^) (£ n r ± ) ] (1-72) 


E 12 = 2 (-l) n+1 w n+1 


jv — —I— 

U n+l U) - + 2 I [(“ 1 ) k<yku k ^) < I n r i) ] -U~ 73 ) 


J 1 - V 

f 21 = + I [(-X) k ^Vl<«' ( 2. nr i’J + KD -1 ( « 


k=n-l 

3=1 
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j-n 

k=0 

F 22 = (-l> n «*J n (€> + I [(-l) kyjk U k (C) (X ? r ± )] + KU o (?) (1-75) 

k=n-l ^ 

j=l 

j=n 

k=2 

G 1 = C-l) n+ V +2 U niWL «) + 'i [(-l) k J\_ 1 (£) (£ n r.)] (1-76) 

k=n+l ^ 

j-1 

j=n 

k=2 

G 2 = (-l) n+2 « n+2 U n+2 (^) +X [(-D^U^U) (£ n r ±)] (1-77) 

k=n+l " ^ 

j-1 

These then are the recursive equations required for the complete 

th 

generalization to a n order system. By employing the above 

equations and replacing in PROGRAM PROJECT cards 100 through 150 

and 300 through 540 with the appropriate programming, the program 

, th 

may be used for any given n order system. 

In like manner by generally defining: 


and: 


P(s) = Q' 2 Ks 2 + apKs + P 2 K 


(1-78) 


J-n 

K=1 


Q ( s ) = s n+2 + 2Ps n+1 + P 2 s n + 2 X (X B r i^ sk + 


k=n 

j=l 


3=n 

k=0 


j=n 

k=2 


s 2 1 £ n r i> s * + 1 <Z ? r i> s * d- 791 

k=n-l - 1 k=n+l ^ 

j=l j=l 

and using Equations (1-49) through (1-52) we may replace in the 
program cards 2860 and 2880 through 2920, thus adapting this 
part of the program to a general n^* 1 order application 
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1.6 COMMENTS 


Throughout this development of the Parameter Plane quadratic 

extension, • tne, c, * s in,tne generalized coefficient form: 
n 

(b k c* 2 +;c£c v + ,<y*p w + e.0 + ,f k P 2 + g ) =0 (1-80)' 


k=0 


i 4 * ■* * * 

have been identically zero. This at first appearance might seem 
to detract from the generalization. The inclusion of this para- 
meter does not however introduce any great difficulty in the 
solution. The change in the development would be to the value 
of P^ and P 2 which would become: 

P x = C 1 + PDj. (1-81) 


P 2 = o 2 + PD 2 


(1-82) 


and the final solution for a which would change to: 


a = 


2E 21 ^ 


*1 - 2 B 21 Q 1 


4B 2 21 


(1-83) 


For this case, new selection rules for acceptable values of or 
would be used, and would be much like those presented for P. 

Though the extension of the Parameter Plane to include the 
cv - P quadratic case makes this tool even more useful, further 
work is still to be done in this field. Not only must the 
equations for the solutions of the Parameter Plane curves for 
such cases as: 


= V 2?2 + =* 


2 

<Y p 


+ d^arp' 




2 . - ~2 
+ V 


g k «p + V + + tjc 


[5] (1-84) 
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and higher ordered combinations of the parameters be developed, 
but more efficient prqgramming techniques must be developed. 

In the use of PROGRAM PROJECT, for instance, as the location of 
the system poles on the a axis of the S-plane move to the left, 
the computational time becomes excessive due to present pro- 
gramming technique and computer speed. 

Another major problem in* further extensions of these tech- 
niques, and indeed even other applications of the curves from the 
proceeding development, will be interpretation. In this case, 
the initial substitution of variables immediate! y allowed inter- 
pretation of the curves sight unseen. Here then, will be most 
likely the one single drawback to further extension, for as the 
parameters a and P are used as representations of other variables 
inricontrol systems, each application will have its own unique 
interpretation . 
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APPENDIX I 


PROGRAM PROJECT is designed to solve the a Quadratic and P 
qua-rtic. The program is divided into two, main sections, the first 
for the computation of the a-p' points and the second for the band- 
width points. 

The first section computes an 80 by 80 matrix of the ot and 
^ points corresponding to set values of £ and w. The computational 
part is followed by two distinct graphing sections* one for lag 
and the other for lead compensation. 

The lag graphing section is set up so that during the plotting 

of the curves each value of or is tested to determine if its value 

-7 

is 10 ^ o; <£ 1..0001. If no points are found within this range 

then a print out is made: 

NO LAG COMPENSATION POSSIBLE 

For the lead, section graphs, & is again tested by the criter- 
ion 1.0001 < ct <, (X-graph scale) (X graph width) . Again if there 
are no values of ® within this region the statement: 

NO LEAD COMPENSATION POSSIBLE 

is printed. In this case however a study of the printed values 
of or must be made to insure that the points are indeed non-existant 

or rather just lie outside the range of the graph. 

* 

The second main section of the program computes the value of 
3 for a given value of ® is determined by the X graph scale. Here 
the plotting routine is set up so as to not plot zero points and 
to stop the curve when either the or or (3 value exceeds the range 
of the graph. 
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PROGRAM PROJECT 

C THIS PROGRAM COMPUTES THE VALUES OF BETA ( POLE LOCATION) AND ALFA (POLE 
C -ZERO RATIO) BASED ON PARAMETER PLANE TECHNIQUES * THE COMPUTED 
C ALFA AND BETA VALUES ARE THOSE REQUIRED TO PLACE THE ROOTS OF ANY 

C THIRD ORDER SYSTEM » TYPE 0,1»2 OR 3» AT A DESIRED ZETA AND OMEGA 

C LOCATION WHILE MAINTAINING A CONSTANT VELOCITY COEFFICIENT* AFTER 
C COMPUTING THE VALUES IT WILL PLOT THE PARAMETER PLANE CONSTANT ZETA 

C CORVES FROM 0*0 TO 0<>9 AND THE CONSTANT OMEGA CURVES FOR EVEY ONE 

C TENTH OF THE VALUE OF THE MAXIMUM VALUE OF OMEGA USED* A 9 BY 15 

C . INCH GRAPH IS OUTPUT BY THE ROUTINE* THIS IS DONE 
C ON TWO SEPARATE GRAPHS » ONE FOR POSSIBLE LAG COMPENSATION AND ONE 

C FOR LEAD COMPENSATION. THE PROGRAM THEN HAS THE ADDITIONAL FEATURE 
C OF COMPUTING AND PLOTTING THE CONSTANT BANDWIDTH CURVES* 

C THE FOLLOWING FEATURES ARE AVAILABLE WITH THE PROPER USE OF THE DATA 
C • CARDS* 

C 1* THE ALFA-BETA COMPUTATIONS MAY OR MAY NOT BE DONE* 

C 2* LAG COMPENSATION MAY OR MAY NOT BE PLOTTED* 

C- 3 o LEAD COMPENSATION 'MAY OR MAY NOT BE PLOTTED* 

XN) C 4* BANDWIDTH COMPUTATIONS MAY OR MAY NOT BE COMPLETED* 

5* BANDWIDTH CURVES MAY OR MAY NOT BE PLOTTED* (AVAILABLE ONLY IP 
• - THE BANDWIDTH COMPUTATIONS HAVE BEEN MADE) . 

THE FOLLOWING DATA CARDS ARE REQUIRED* 

**CARD ONE “ A> B »C»G - TEN COLUMNS PER NUMBER IN FLOAT 
y THESE ARE THE LOCATIONS OF THE UNCOMPENSATED POLES AND THE 
UNCOMPENSATED ROOT LOCUS GAIN. 

**CARD TWO - WFIN - TEN COLUMNS IN FLOATING POINT. 

THIS IS THE MAXIMUM VALUE OF OMEGA TO BE USED* 

**CARD THREE - IABCMP - COLUMN ONE IN FIXED POINT 

0 - THE ALFA-BETA COMPUTATIONS WILL BE DONE 

1 - THE ALFA-BET A COMPUTATIONS WILL NOT BE DONE 




C 

c 

c 

c* 

C 

C 

c* 

c 

c* 

c 

c 

c 

c 

c 


IF IABCMP-1 CARDS FOUR THROUGH THIRTEEN ARE OMITTED 

C***CARD FOUR ~ ILGPLT - COLUMN ONE IN FIXED POINT 
C 0 - THE LAG ZONE CURVES WILL BE PLOTTED 


000000 



C 1 - THE LAG ZONE CURVES WILL NOT BE PLOTTED 

C ****** 

C -IF ILGPLT=1 THE NEXT FOUR CARDS ARE OMITTED 
C • ’ . ' ' ■ ****** 

C***CARD FIVE ~IT( 1 )~IT( 6 ) - COLUMNS 1-48 IN ALFANUMERIC CHARACTERS 
C this is THE FIRST LINE OF THE LAG GRAPH TITLE 

C***CARD SIX - ITC7) — IT (12) - COLUMNS 1-48 IN ALFANUMERIC CHARACTERS 
C, THIS -IS THE SECOND LINE OF THE LAG GRAPH TITLEo 

C***CARD SEVEN - LBL ( 1 1 ) -L 8 L { 20 J - FOUR COLUMNS PER LABEL (TEN LABELS IN 
C CONSECUTIVE COLUMNS) IN ALFANUMERIC CHARACTERS* 

C THESE ARE THE LABELS .TO BE PUT ON THE CONSTANT OMEGA CURVES* TO 
C DETERMINE WHICH VALUES WILL BE PLOTTED 9 DIVIDE WFIN BY 10 o THIS 

C VALUE AND INTEGER MULTIPLES OF IT TO 10 WILL BE PLOTTEDo 

C***CARD EIGHT - XLGZsYLGZ - TEN COLUMNS PER NUMBER IN EXPONENTIAL OR ' 

C FLOATING POINT 0 . . 

C THESE ARE THE X AND Y SCALES FOR THE LAG GRAPH*, ONLY ONE SIGNI-. 

C- FICANT NUMBER IS. TO BE USED, s 

C***CARD NINE.- ILDPLT - COLUMN ONE IN FIXED POINT 
C 0 - THE LEAD CURVES WILL BE PLOTTED 

C 1 t THE LEAD CURVES WILL NOT BE PLOTTED 

C . ****** 

C IF I LDPLT=1 THE NEXT FOUR CARDS ARE OMITTED 
C • ****** , 

C***CARD TEN - THE SAME AS CARD FIVE EXCEPT FOR mt LfcAD GRAPH 
C***CARD ELEVEN - THE SAME AS CARD SIX EXCEPT FOR THE LEAD GRAPH 

C***CARD TWELVE - THE SAME AS CARD EIGHT EXCEPT FOR THE LEAD GRAPH 

C***CARD THIRTEEN - A***DUPL ICATE*** • OF CARD SEVEN 
C***CARD FOURTEEN - I 8 WCMP - COLUMN ONE FIXED POINT ■ 

C 0 - BANDWIDTH COMPUTATIONS AND GRAPHING WILL NOT BE DONE® 

C 1 - BANDWIDTH COMPUTATIONS WILL BE DONE 

C ****** 

C IF I8WCMP=0 THE REMAINING CARDS ARE OMITTE 
C ****** ■ , 

C***CARD FIFTEEN - BWXjBWY - THE SAME AS CARD EIGHT EXCEPT FOR THE 
C BANDWIDTH CURVES* 

C BWY IS ALSO USED TO DETERMINE WHICH VALUES OF ALFA WILL BE USED IN 



"THE BANDWIDTH COMPUTATIONS,, 

***CARD SIXTEEN - WEND - TEN COLUMNS IN FLOATING POINT 

THIS IS THE MAXIMUM VALUE OF OMEGA FOR WHICH THE BANDWIDTH 
COMPUTATIONS WILL BE DONE 

*#*CARD SEVENTEEN - IBWPLT - COLUMN ONE IN FIXED POINT 

0 - THE BANDWIDTH CURVES WILL BE PLOTTED 

1 - THE BANDWIDTH CURVES WILL NOT BE PLOTTED 

IF IBWPLT=1 THE REMAINING CARDS ARE OMITTED 

***CARD EIGHTEEN ~ THE SAME AS CARD FIVE EXCEPT FOR THE BANDWIDTH CURVES 
**#CARD NINETEEN “ THE SAME AS CARD SIX EXCEPT FOR THE BANDWIDTH CURVES 
*##CARD TWENTY ~ BANDWIDTH CURVE LABELS 

TO DETERMINE WHICH CURVES WILL BE PLOTTED? DIVIDE WEND BY 10 * 

THE PROGRAM PLOTS THIS CURVE AND INTEGER MULTIPLES OF IT UP TO 10 


IT IS RECOMMENDED THAT FOR THE INITIAL RUN THE FOLLOWING DATA CARDS 
BE USEDo 

CARDS 1 9 2 s> 3 ( I ABCMP = 0 ) s>A{ ILGPLT = 1 ) ?9( ILDPLT = 1 ) 9 1A( IBWCMP=0) 

THESE DATA CARDS WILL ALLOW ONLY THE ALFA-BETA COMPUTATIONS IU bt 

COMPLETEDo A PRINT OUT OF THE VALUES WILL BE OUTPUT WHICH WILL ALLOW 
YOU TO CHOOSE THE PROPER CURVES AND SCALESo CAREFUL SELECTION 
OF CURVE SCALES IS IMPORTANT ? FOR THE PROGRAM WILL NOT ALLOW POINTS 
OUTSIDE THE AXIX LIMITS TO BE PLOTTEDo 


DIMENSION AFIN(80?80) 9 BFIN<80 9 80> 9 XAZ(80I *YBZ(80) jXAW(80) » 000010 

1 YBW{ 80 ) » IT (12) sLBL{ 20 ) 9 BC0FI (5) *R00TR(4) s ROOT I (4) 9 ACOFI (3) 9 000020 

2 U( 101 9 AROOTl (A) 9 ACOFR<3) ?BCOFR(5) 9 WLAB ( 80 ) sZLAB { 80 } 9 AROOTR (A) 000030 

COMMON BCOFR 9 BCOFI jROOTR « ROOT I 9 BFI NAL » I FLAG 9 AFI N 9 BF I N 0000A0 

9999 PRINT 1A0 000050 

1A0 FORMAT (1H1) 000060 

DO 60 JK = 1 » 6A00 000070 

AFlNtJK) =0o0 000080 

60 BFINCJK) = 0*0 000090 



READ IjAjBoCsG 

1 FORMAT (4F 10 <,0 5 
PROD = A*B*C 
SUM = A+B+C 

o.MPRD = A*B 4- A*C + B*C 
PRDGN = PROD + G 
ZETA = OcO 
READ 2 oWFIN 

2 FORMAT ( F10 o 0 5 
READ 9 1 IABCMP 

9 FORMAT (II) 

I F < I ABCMP-1 5 23 9 24 9 24 
23 STP = WFIN/80 <s 
DO 12 L = 1*80 
LJ « 80* { L~1 ) 

W « STP 
30 U<1)=-1* 

U { 2 5 =0 o 
U( 3 ) =1 o 
DO 10 N = 2 s 6 

U<N*2) =2o*ZETA*U(N+l)-U(N) 
J DO 11 J=1 ®80 

^ LJ = LJ. + 1 

V W2=W*W 

W3=W2*W 
W4=W2*W2 
W5»W2*W3 
CONN = G*W2 
CON « -2**G*W 
CONI = 2 o*W4 
CON2 = -2 o*SUM*W3 
C0N3 = 2o*SMPRD*W2 
CON4 = -2,*PROD*W 
CON 5 « SUM*W2 ’ 

C0N6 =-~SMPRD*W 
CON7 - SUM*W4 


000100 

000110 

000120 

000130 

000140 

000150 

000160 

000170 

000180 

000182 

000185 

000188 

000190 

000200 

000210 

000220 

000230 

000240 

000250 

000260 

000270 

000280 

000290 

000300 

000310 

000320 

000330 

000340 

000350 

000360 

000370 

000380 

000390 

000400 

000410 

000420 






C0N8 = -*SMPRD#W3 
CON 9 * PROD*W2 
B21 = C0NN*U<3) 

B22 = C0NN*U(4) 

D1 « CON*U‘<2) 

D2 = CON^U ( 3 ) 

•■£11 = CONl^U ( 5 ) + CON2^U ( 4 ) + CON3*U(3) + C0N4*UJZ) 

E12 = CONl*U<6) + CON2*U( 5 ) + CON3*U(4) + CON4*U(3) 

F21 = -W3#U<-4) + CON5*U(3) + CON6*U(2) + PRDGN#U(1) 

-F22 = “W3*U(5) + CON5*U(4) + CON6*U(3) + PRDGN*U(2) 

G1 = -W5*Ut6) + CON7*U(5) + CON8*U(4) + CON9*Ut3) 

62 = -W5*U<7) + CON7*U<6) + CON8*U(5> + CON9*U(4) 

COFl = B21*F22*(2o*F2l*B22-F22*B21)-'F21*(F21*B22*B22+D2#D2*B2l} 
COF2 = Ell* ( 2«*B22*( F22*B21-F21*B22 > ~D2*D2*B21 )+2 o*E12*B21* ( F21*B2 
12-F22*B21) ' ‘ • * ' 

COF3 •«= B21*<“B21*lE12*E12+2**F22*G2)-D2#D2*61+2o»822*(Ell»E12+F21* 
1 G2+G1*F22) )-B22*B22*<Ell*Ell+2**F21*Gl > 
COF4=2**G2*B21*(Ell*B22-E12#B21)+2o*B22*Gl*( E12#B23-F1 1*B22) 

• COF5= -(G2*B21-Gl#B??)*fG2*B21-Gl#B22i 


DO 30 I =1 » 5 
50 BCOFI ( I) = OoO 
BCOFR(l) - loO 
BCOFR ( 2 ) = COF2/COF1 
8COFR { 3 ) = COF3/COF1 
.BCOFR ( 4 ) ~ COF4/COF1 
BCOFR ( 5 ) “ COF5/COF1 
CALL ABETART 
IFLAG = 0 
CALL SORT 

IF (IFLAG-1)300»11»11 
500 BFIM(LJ) = BFINAL 

Q1 = BFIM<LJ)*(E11+BFIN(LJ)#F21)+G1 


ACOFR ( 1 ) "lo 0 
ACOFR ( 2 ) = Oo 0 
ACOFR < 3 ) = Q1/B21 
ALFASQ « ABSF ( ACOFR (35) 


000430 

000440 

000450 

000460 

000470 

000480 

000490 

000500 

000510 

000520 

000530 

000540 

000550 

000560 

300570 

300580 

300590 

000600 

000610 

000620 

000630 

000640, 

000650 

000660 

000670 

000680 

000690 

000700 

000710 

000720 

000730 

000740 

000750 

000760 

000770 

000780 



AFIN(LJ) = SQRTF ( ALFASQ ) 

11 V/ = W+STP 

12 2ETA = ZETA + o 0125 
LBL(l) = 4HZ=o0 
LBL ( 2 ) = 4HZ"o 1 
LBL<3) = 4HZ= o 2 


LBL (4) = 4HZ=o3 
LBL (5) «' 4HZ= o 4 
LBL ( 6 ) = 4HZ~ o 5 
LBL (7) = 4HZ= o 6 
LBL (8) = 4HZ=o7 
LBL (9) = 4HZ=c8 
LBL ( 10 ) « 4HZ = o9 
READ 7 j ILGPLT 

7 FORMAT (II) 

IF C ILGPLT-1)8»67»67 

8 READ - 3 s' ( IT C I) » 1 = 1 5 12 > 

3 FORMAT (6A8) 

READ 6 i ( LBL ( I ) pl-ll»20>- 
6 FORMAT ( 10A4 ) 


READ 4s' XLGZ jYLGZ 
4 FORMAT (2E10a0) 

XLGLM = 9«*XLGZ 
YLGLM «' 15o*YLGZ 
MODE = 1 
IL = 0 

DO 62 K=ls8Uja 
LL = 1 

KJ = ( K~1 }*80 
DO 61 J = 1 3 80 " 

KJ ~ KJ+1 

IF (AFINtKJ 5-o 0000001 ) 61, 61 10 *6110 
6110 IF(AFIN(KJ)“lo0001 ) 6113*61 9 61 

6113 IF(AFINtKJ) - XLGLM) 6114»61 961 

CARDS 1120 - 1130 ARE MISSING 

6114 XAZ(LL) = AFIN(KJ) 


000790 

000800 

000810 

000820 

000830 

000840 

000850 

000860 

000870 

000880 

000890 

000900 

000910 

000920 

000930 

000940 

000950 

000960 

000970 

000980 

000990 

•001000 

001010 

001020 

001030 

001040 

001050 

001060 

001070 

001080 

001090 

001095 

001100 

001110 


001140 



001150 


I F { BF I N < KJ } - YLGLM )61 12*61*61 


CARDS 1160 - 1170 ARE MISSING 


6112 

YBZ(LL) = • BF I N ( KJ ) 

001180 


LL * LL + 1 

001190 

61 

CONTINUE 

001200 


JJ s LL - 1 

001210 


IL- = IL + 1 

001220 


I F { JJ-l )62»62»6116 

001230 

6116 

LAL = LBL< IL) 

001240 


CALL DRAW( JJ*XAZ»YBZ*MODE»0*LAL* IT »XLGZ*YLGZ *0*0*0* 0*9* 15*0* LAST) 

001250 

6111 

MODE « 2 

001260 

62 

CONTINUE 

001270 


IF<MODE~1)65*65»6120 

001280 

6120 

DO 66 K =8 9 80 s 8 

001290 


LL * 1 

001300 


DO 63 J=l*80 

001310 


JK = ( J~1 )* 80 + K 

001320 


IF<AFIN( JK)-«0000001 ) 63 ,6127 * 6127 

001325 

6127 

IF(AFIN( JK)~lo 0001 ) 6123 *63*63 

001330 

6123 

I F ( AFI N ( JK ) - XLGLM)6124*63*63 

001340 

C 

CARDS 1350 - 1360 ARE MISSING 


6124 

XAW(LL) = AFIN(JK) 

001370 


IFIBFIN(JK) - YLGLM ) 6122 * 63 * t>3 

001380 

C 

CARDS 1390 ~ 1400 ARE MISSING 


a 6122 

YBW(LL) = BFIN(JK)' 

001410 

K-) 

LL = LL + 1 

001420 

(S' 63 

CONTINUE 

001430 


JJ = LL - 1 

001440 

IL = IL + 1 

001450 


IFIJJ- 1)6121*6121*612* 

001460 

6121 

I F < K-80 >66*6125*6125 

001470 

6125 

MODE * 3 

001480 


LAL = 4H 

001490 


JJ = 2 

001500 


XAW ( 1 ) = XLGLM 

001510 


XAW ( 2 ) = XLGLM 

001520 



YBW(l) * 0*0 
YBW ( 2 ) = YLGZ 
GO TO 2000 
6126 LAL = LBL(IL) 

2000 CALL DRAW< JJ*XAW*YBW»MODE»0»LAUIT*XLGZ»YLGZ*0»0 9 0*0*9*15*0*LA$TI 
MODE « 2 
IF{K~7" )66?64?64 

64 MODE = 3 

66 CONTINUE 
GO TO 67 

65 PRINT 130 

.30 FORMAT ( IX ? 33H NO LAG COMPENSATION IS POSSIBLE •>///) 

67 READ 20? ILDPLT 

20 FORMAT (II) 

IF{ ILDPLT“l)68?10uu?iuuu 

68 READ 5? < IT ( I) ? 1 = 1 ? 12 ) 

5 FORMAT (GAS) 

READ 21? XLDZ?YLDZ 

21 FORMAT (2E10„0) 

READ 22 9 ( LBL C I ) }I=ll»20 

22 FORMAT (10A4) 

XLDLM = 9 U *XLDZ 
YLDLM = 15o*YLDZ 
I L a 0 

MODE a 1 

DO 72 !< = 1 ? 80 ? 8 
KJ = ( K-l > *80 
KK = 1 

DO 71 J = 1 9 80 
KJ = KJ + 1 

IF(AFIN(KJ)“lo0001 ) 71? 71 11 ?7111 
7111 IF ( AFIN ( KJ ) - XLDLM)7ll7»71»7l 
C CARDS 1830 - 1840 ARE MISS'ING 

7117 XAZ(KK) a AFIN(KJ) 

~F(BFIN( KJ ) - YLDLM)7118?71>71 
ARDS 1870 ~ 1880 ARE MISSING 


001530 

001540 

001550 

001560 

001570 

001580 

001590 

001600 

001610 

001620 

001630 

001640 

001650 

001660 

001670 

001680 

001690 

001700 

001710 

001713 

001716 

001720 

001730 

001740 

001750 

001760 

001770 

001780 

001790 

001800 

001810 

001820 

001850 

001860 



/z° 


7118 

YBZ(KK) - BFIN(KJ) 

001890 


KK « KK + 1 

001900 

71 

CONTINUE 

001910 


MM = KK~1 

001920 


IL = IL + 1 

001930 


IFCMM-D72 *72 •7119V 

001940 

7119 

LAL = LBL ( IL) 

001950 


CALL 0RAW(KM»XA7.YR7 tMODE»OsLALj IT sXLDZs YLDZ »0 »0 * 0 • 0 >9> 15 »0 » LAST ) 

001960 

7110 

MODE * 2 

001970 

72 

CONTINUE 

001980 


IF (MODE-1 >70*70 • 78 

001990 

”8 

DO 76 K=8 p 80 <> 8 

002000 


KK * 1 

002010 


DO 73 J « 1 9 80 

002020 


JK = {J-l)#80 + 

002030 


I F{ API N(JK)-loOOOl) 73*7121 »7121 

002040 

7121 

IFtAFIN(JK) - XLDLM)7127*73»73 

002050 


CARDS 2060 - 2070 ARE MISSING 


71 71 

XAV/(KK) , = AFIN( JK) 

002080 


I F ( BF I N ( JK ) - YLDLM ) 7128 9 73 9 73 

002090 


CARD'S -2100 ~ 2110 ARE MISSING 

.... 

7128 

YBW(KK) ~ BFIN(JK) 

002-120 


KK = KK + 3 

002130 

73 

CONTINUE 

002140 


MM = KK-1 

002150 


IL «■ IL + ] 

002160 


I F ( MM- 1 ) 7 1 2 u » a’20»712<> 

002170 

7120 

IF(K-80)76*7122>7122 

002180 

7122 

MODE « 3 

002190 


LAL - AH 

002200 


MM * 2 

002210 


XAW(l) - ' XLDLM 

002220 


XAW i 2 ) « XLDLM 

002230 


YBW( 1) ' a OoO 

002240 


YBW ( 2 ) ■ * 'YLDZ 

002250 


GO TO 2001 

002260 



7129 

2001 


75 

76 

70 

131 

1000 


81 ■ 

32 
100 
101 
102 
103 
111 
i. JL C 



219 


CALL DRAW ( MmIxAW * YBW 9 MODE > 0 s LAL » 1 T *XLDZ * YLDZ » 0 * 0 * 0 » 0 *9 * 1 5 * 0 * LAST ) 
MODE = 2 

IF(K-72>76975975 
MODE = 3 
CONTINUE 
GO TO 1000 
PP I NT 131 

FORMAT (IX « 34H NO LEAD COMPENSATION IS POSSIBLE t///) 

CONTINUE 

ZLAB(l) = 0 o 0 

DO 81 1=2910 

ZLAB ( I ) ' = ZLABC 1-1 ) + 

WLABC8) = 8c*STP 

DO 82 N = 16»80<>8 

WLAB(N) = WLAB(N-85 + WLAB ( 8 ) 

PRINT 100 
FORMAT ( 1H1 ) 

PRINT 101 • 

FORMAT! 2Xs 2 OH THE ALLA VALUta A^ti 9 / ' > 

PRINT' 102 9 (ZLAB ( I ) I = X s> 1 0 ) 

FORMAT UX 96 H ZETA slOFllo 6 > 

PRINT 103? {WLAB( J ) 9 ( A/IN ( J 9 1 ) 9 1=1 973 9 8 ) J = 8 9 8 U 9» ) 

FORMAT (/9lX9F6c29l0Ello5) 

PRINT 111 

FORMAT (///// 2X *20 H THE BE I A VALULi. arc 
PRINT 1129 , (ZLABtl) 1=1*10) 

FORMAT (1X*6H ZETA '*10Fll*6l ■ 

PRINT- 113® ( WLAB ( J ) * ( BFIN ( J 9 1 ) 9 1 = 1 9739 8 ) J= 8 y 80 98 

FORMAT (/*lX*F6o'2*10Ell*5) 

PRINT 114 
FORMAT ( 1 H 1 ) 

READ 2189 IBWCMP 
FORMAT (ID 

IF ( IBWCMP-1 ) 1002 s-219 9 219 
READ 217 ? BWX 9 BWY 


002270 

002280 

002290 

002300 

002310 

002320 

002330 

002340 

002350 

002360 

002370 

002380 

002390 

002400 

002410 

002420 

002430 

002440 

002450 

002460 

002470 

002480 

002490 

002500 

002510 

002520 

002530 

002540 

002550 

002560 

002570 

002580 

002590 

002600 

002610 

002620 



217 FORMAT (2E10o0) 002630 

READ 223 ? WEND 002633 

223 FORMAT (FlOoO) 002636 

YBWLM «' 15»*BWY 002640 

XLM * 9#*BV/X 002650 

AM2*«5 002660 

STEP = WEND/20* 002670 

W - STEP 002680 

ALFASP « XLM/20* 002690 

XAWU) = ALFASP 002700 

DO 200 K=2« 20 002710 

200 XAV/ ( K } = XAW(K-l) + XAW(l) 002720 

DO 203 N=1 s 20 002730 

ALFA = ALFASP 002740 

DO 210 M=l»20 002750 

XA2 ( M) « OoO 002760 

210 Y82 ( M ) - OoO 002770 

DO 207 1=1 » 20 002780 

YBW(N) ** W 002790 

W2**W#W 002800 

W3“W*W2 002810 

W*=W2*W2 002820 

W5*W2#W3 002830 

A6R=,-W2-v 002840 

AKR * 6 002850 

ADR=2o*W4~ 2*#W2#SMPRD 002060 



All = 2**G*W 

AER=“V/2*S0M + PROD 

AFR=l'/4*SUM - W2*PR0D 

ADI -~2o*W3XSUM + 2o#W*PR0D 

AEI « -W3 + W*SMPRD 

AFI = W5 ~ W3*SMPRD 

PI « ADR*ADR+ADI*ABI 

Q1 = 2 ©*ADR#AER + 2©#AEI*ADI 

R1 = 2o*ADR*AFR + AERfcAER + AEI*AEI + 2*#ADI*AFI 
''1 = 2<j#AER*AFR + 2o*AEI*AFI 


002870 

002880 

002090 

002900 

002910. 

002920 

002930 

002940 

002950 

002960 




W1 = AFR*AFR + AFI^AFI 
A2 = ALFA*ALFA 
A 4 = A2*A2 
R2 = A2*AII#AI 

WW2 = A4*AGR*AGR +. 2 o *A2*AKR#AGR + 
DO 51 M = 1,5 
51 BCOFI(M)„= OoO 
BCOFR { 1 ) = loO 
BCOFR ( 2 ) = Ql/Pl 
BCOFR ( 3 ) « <AM2*R1~R2)/(AM2*P1) 
BCOFR (4) = Vl/Pl 
BCOFR ( 5 ) = ( AM 2 *W 1 -WW 2 ) / ( AM2 #P 1 ) 
CALL' ABETART 
IFLAG = 0 
CALL SORT 

IF(IFLAG-1)900 9 206 5 206 

206 BFINCN* I) ■ * OoO 
GO TO 207 

900 BFINtNsU « BFINAL 

207 ALFA « ALFA + ALFASP 

203 W *.'W + STEP 
READ 216 9 IBWPLT 

216 FORMAT (III 

IF ( IBWPLT-1 ) 214 s 1001 9 1001 
214 MODE « -1 

READ 202 ? ( r I T { K ) 9K = lol2 ) 

202 FORMAT (6A8) 

READ- 201® (LBL(N) 9N=2s>20«.2) 

201 FORMAT <10A4> 

DO 211 N = 2 9 20 9 2 
KK- a 1 

IF(N-20)204*205*205 
205 MODE a 3 

204 CONTINUE 

DO 212 1 = 1 9 20 

IF(BFIN(N*I ) o 00000 1 ) 212 9 2 12 9 209 


002970 

002980 

002990 

003000 

AKR*AKR 003010 

003020 

003030 

003040 

003050 

003060 

003070 

003080 

003090 

003100 

003110 

003120 

003130 

003140 

003150 

003160 

003170 

003180 

003190 

003200 

003210 

003220 

003230 

003240 

003250 

003260 

003270 

003280 

003290 

003300 

003310 

003320 



003330 

003340 

003350 


209 I F ( BFIN (N $ I ) ~ YBWLM ) 905 > 906 » 906 

905 YBZ(KK) = BFIN{N>I ) 

XAZ(KK) = XAW(i-) 

CARD 3360 IS MISSSIN6 
KK. = K.K + 1 
212 CONTINUE 

906 JJ = KK -s 1 

IFC JJ-1 )221 #221*^u 
221 IF(N-20)211*225»225. 

225 I F { MODE-1 )224j224s>215 
15 MODE =' -3 
LAL ’ - 4H 


220 

2002 

222 

208 

211 


XAZ ( 1 5 = 0 o 0 
XAZ (2) = 0*0 
YBZ(l) a OoO 
YBZ < 2 ) = BWY 
JJ = 2 
GO TO' 2002 

LAL = LBL(N) . A 

CALL DRAW ( JJ*XA2 » YBZ »M0DE*0»LAL* IT»BWX»8V/Y*O*0*0 >0>9 # 15»0 *LAST ) 

IF ( N-20 }208o2Tl$211 

MODE = 2 

CONTINUE 

GO TO 1001 


226 FORMAT < 1H1 » IX 9 76HTHERE ARE -NO POSSIBLE BANDWIDTH CURVES FOR THE F 
' 1INAL VALUE OF OMEGA STATED 
GO TO 1002 

1001 PRINT 120 

120 FORMAT ( 1H1 ® 20X » 38H VALUES OF BETA FOR CONSTANT BANDWIDTH 
.•PRINT 122s (YBW(K) 9 K=*2s20,2) 

122 FORMAT ( / « 9X »10F11 06 > 

PRINT 121 ? <XAW(K) » ( BFINC J»K) » J«2920>2) 9K«ls20 J 

121 FORMAT (/»1XjF6o2s2Xj10E11o6) 

1002 CONTINUE 
CO TO 9995 


003370 

003380 

003390' 

003400 

003410 

003415 

003420 

003430 

003440 

003450 

003460 

003470 

003480 

003490 

003500 

003510 

003520 

003530 

003540 

003541 

003542 

003543 

003544 

003545 

003550 

003560 

003570 

003580 

003590 

003600 

003610 

003615 



END 



1 jCONV ( 50 ) 

DIMENSION AF IN (80 9 80) 9 3FIN{80 5 80) 

COMMON A 9 Y I MAG 5 U 9 V 9 DUMMY 1 9 DUMM Y2 9 AF I N 9 3 FI N 
M = 4 

= 10o0 

.-25 

ER=0 

*F(M) 54o54952 
54 I ER=1 
52 NP3=N+3 

100 B( 2 )~0o0 
B( 1 ) “0 0 0 
C ( 2 ) =0* 0 
C(l>=0 o 0 
D ( 2 ) =Oo 0 
E ( 2 ) =0o 0 
H( 2 ) =0o.0 

DO 101 J=3 e NP3 

101 H(J)=A(J~2) 

T = 1 oO 

SK*10o0**F 

150 I F ( H ( NP3 5 1 2009 151 ?200 
15 KNP3) = 0<,0 
^(NP3)*0«0 
CONV<NP3)»SK 
NP‘3=NP3-1 

IF(NP3) 15291529150 
152 I ER"1 ’ 

200 IF(NP3“3)205 9 5l920l 
205 I ER=1 

201 P$=0o0 


003620 

00363 

00364 
003650 
003660 
003670 
003680 
003690 
003700 

00003710 

00003720 

00003730 

00003740 

00003750 

00003760 

00003770 

00003700 

00003790 

00003800 

00003810 

00003820 

00003830 

00003840 

00003850 

00003860 

00003870 

00003880 

00003890 

00003900 

00003910 

00003920 

00003930 

00003940 

00003950 

00003960 

00003970 



QS=OoO 

PT=OoO 

QT=OoO 

S=OoO 

REV? lo 0 

SK-10o0**F 

IF(NP 3-4)206920252 03 

206 IER = 1 

202 R-~H ( 4 ) /H ( 3 ) ’ 

• GO TO 500 

203 DO 207 J=3 9 NP3 

• IF(H( J) )204 ? 2079204 
204' S=S+LOGF(ABSF<H( J)‘) ) 

207 CONTINUE 
FPNl=N-i-l 
S=EXPF ( S/FPihx j 
DO 208 J=3 « NP3 

208 H(J)=H(J)/S 

210 I F { ABSF 4 ) /H { 3 ) ) -ABSF(H ( NP3-1 ) /H ( NP3 ) ) )250 9 2529252 

250 T=-T 

( NP3-4) /2 + 3 
DO 251 J=3 9 M 
S«H ( J ) 

JJ=NP3~J+3 
H(J)=H< JJ) 

251 H ( J J ) =S 

252 IF(QS) 25392549253 

253 P=PS 
Q=QS 

GO TO 300 

254 HH2~H( NP3-2 ) 

IF(HH2) 256 p255 9 256 

255 0=lo0 
P=-2o0 

GO TO 257 

256 Q=H ( NP3 ) /HH2 


00003980 
00003990 
00004000 
00004010 
00004020 
00004030 
00004040 
00004050 
00004060 
00004070 
00004080 
00004090 
00004100 
00004110 
00004120 
00004130 
00004140 
00004150 
00004160 
00004170 
00004180 
00004190 
00004200 
00004210 
00004220 
00004230 
00004240 
00004250 
00004 260 
00004270 
00004280 
00004290 
00004300 
00004310 
00004320 
00004330 



P= (H<NP3-l)*-Q*H(NP3-3) )/HH2 
IF ( NP3-5 ! 258*5509258 
R=0«0 

DO 490 1=1 »L 

DO 351 J = 3 s NP3 

B ( J ) =H ( J ) -P*B ( J-l ) -Q*B ( J~2 ) 

C ( J ) =B ( J } ~P*C { J~ 1 ) ~Q*C ( J-2 ) 
IF{H(NP3“1 J ) 352 s400 s 352 
IF ( B ( NP3-1 ) ) 3 53 s 400s 353 
AVHB1=ABSF(H(NP3“1 )/B(NP3-l ) ) 
IFCAVH81~SK)450s354s354 
B ( NP3 ) =H ( NP3 ) -Q*B ( NP3-2 ) 

IF ( B ( NP3 ) ) 401 s 5 50 s 401 
AVHB2=ABSF { H ( NP3 ) /8 ( NP3 ) ' 
IF(SK-AVHB2)550s450s450 
DO 451 J = 3 s NP3 
D( J)=H( J)+R*D< J-l) 

E( J)=D( J)+R*E( J-l) 

IF ( D ( NP3 ) ) 452 s 500 s 452 
AVHD3=ABSF (H(NP3)/D(NP3) ) 

IF ( SK— AVHD3 ) 500 s453s453 
CC2=C ( NP3-2 ) 

CC3=C ( NP3-3 ) 
aNP3-l)=~P*CC2-G*CC3 
CC1=C ( NP3-1 ) 

S=CC2*CC2-CC1*CC3 
IF ( S ) 455 s 454 9 455 
P=P-2o0 • 

Q=Q* (Q+l o 0 ) 

GO TO 456 

P= P+ ( B ( NP3-1 1 *CC2~B ( NP3 ) *CC3 ) /S 
Q=Q+{-8(NP3-l )*CC1+B(NP3)*CC2)/S 
I F ( E ( NP3-1 ) ) 4 5 8 9 4 5 7 s 4 5 8 
R=R-loO 
GO TO 490 

R=R~D(NP3)/E?NP3~1 j 




00004340 

00004350 

00004360 

00004370 

00004380 

00004390 

00004400 

00004410 

00004420 

00004430 

00004440 

00004450 

00004460 

00004470 

00004480 

00004490 

00004500 

00004510 

00004520 

00004530 

00004540 

00004550 

00004560 

00004570 

00004580 

00004590 

C 3004600- 

: 0004610 

00004620 

00004630 

00004640 

00004650 

00004660 

00004670 

00004680 

00004690' 



490 

CONTINUE 

00004700 


PS=PT 

00004710 


QS=GT 

00004720 


PT =P 

00004730 


QT^Q 

00004740 


I F ( REV) 491 s 492 j 492 

00004750 

491 

SK=SK/10. 0 

00004760 

492 

REV=-REV 

00004770 


GO TO 250 

00004780 

500 

I F ( T ) 50 1 s 502 s> 502 

00004790 

501 

R=1 «>0/R 

00004000 

502 

NP=NP3-3 

00004010 


U(NP)=R ' 

00004020 


V( NP)“0o0 

00004830 


CONV(NP ) =SK v. 

00004840 


NP3=NP3-1 

00004050 


DO 503 J=3 9 NP3 

00004860 

503 

H( J)=D( J) 

00,0;3-4870 


IF(NP3~3)300$i5l5300 

00004080 

550 

I F { T ) 551 ?552 « 552 

00004890 

551 

P=P/Q 

00004900 


Q=1 oO/Q 

00004910 

552 

PP2=P/2oO 

00004920 


QMPSQ~Q~PP2*PP2 

00004930 

560 

IF(QNPSQ) 554 >5545 553 

00004940 

553 

NP=NP3-3 ' 

00004950 


U{ NP }=~PP2 

00004960 


U( NP-l ) =~PP2 

00004970 


S--SQRTF IQMPSQ) 

00004980 

fj 

V ( NP ) =S ' 

00004990 

[(• i 

V( NP-l ) *="S 

00005000 

V\\ 

GO TO 561 

00005010 

554 

S«SQRTF(~QMPSQ) 

00005020 


NP=NP3-3 

00005030 


IF(P)555 ! >556 s 556 

00005040 

555 

U( NP )="*PP2+S 

00005050 



GO TO 357 
556 U(NP)=-PP2-S 
357 U ( NP-1 ) “Q/U ( NP ) 
■ V( NP ) =0 oO 
V ( NP~1 ) “0 oO 
561 CONV ( NP ) = SK. 
CONV ( NP~1 ) =SK 
NP3=NP3~2 
DO 558 J = 3 o NP3 
558 H(J)=B(J) 

GO TO 200 
51* RETURN 
■ END 


SUBROUTINE SORT' . 

DIMENSION RR(4) *RI (4) 9 REAL ( 5 ) * Y I MAG ( 5 ) 
DIMENSION AFINt 80' >80 ) >BFIN(80o80 ) 

COMMON REAL * Y I MAG 9 RR ?R I # SF « I FLAG 9 AFIN 9 BF IN 
B = OoO 
DO 800 1*1*4 

IF ( ABSF (RI( I ) )“loE~7)801»800s>800 
801 B = MAX IF ( B >RR ( I) ) 

800 CONTINUE 

I F { B ) S02c>802?803 
80 2 PRINT 8 04- 

804 FORMAT ( 34H THERE ARE NO Pubilivt klal. kwio 
IFLAG = 1 
RETURN 
803 BF = B 
RETURN 
END 
END 



N 


00005060 

00005070 

00005080 

00005090 

00005100 

00005110 

00005120 

00005130 

00005140 

00005150 

00005160 

00005170 

00005180 

00519 

00520 
005210 
005220 
005230 
005240 
005250 
005260 
005270 
005280 
005290 
005300 
005310 
005320 
005330 
005340 
005350 
005360 
005370 
005300 




Fig. 1-1. Double Lead Compensation of Plant 
with G (s) = 
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Fig. 1-5. Double Lead Compensation of Plant 
with G (s) 
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TRANSIENT RESPONSE OF NONLINEAR SYSTEMS 
2 . 1 INTRODUCTION 

When a system has one nonlinear element that is single val- 
ued and non-frequency dependent, analysis of the system is con- 
veniently accomplished using the parameter plane methods. The 
nonlinear element is represented by a describing function, which 
is a function of signal amplitude only. The describing function 
is designated as one of the parameters, ct or /3. This designation 
removes the nonlinear parameter from the functions that determine 
the parameter plane curves so that these may be plotted on the 
a-j3 plane. The M-point is located on the ce-jS plane in the usual 
way, but for the case of one nonlinear element one coordinate of 
the M-point is the numerical value of the describing function of 
the nonlinear parameter. I^or linear systems the M-point is sta- 
tionary on the a-jS plane, but for a nonlinear system the M-point 
moves because the numerical value of the describing function is 
a function of signal amplitude. For a system with one single 
valued nonlinearity, N, where N is designated as j8, the locus 
followed by the M-point is a straight line parallel to the /3-axis. 
This locus of M-point motion can be said to start at the value of 
|3 corresponding to very small (zero) signal amplitude into the 
nonlinear element. The displacement of the M-point along this 
locus is determined by the way in which jS varies as a function 


.of signal amplitude, and this is determined by using the 


Constant -£ and constant -co 
-co curves. 



or constant -cr and constant 


the describing function of the nonlinear element. 

Previous work has shown how to predict limit cycles using 
M-point locus on the parameter plane. If this locus crosses the 

i 

stability boundary (C=0 curve or <?=0 curve) the intersection of 
these curves defines the frequency of the limit cycle. If an 
amplitude scale can be determined for the location of the M-point 
on the M-locus;, then this scale is used to define the amplitude of 
the limit cycle. 

The concept of a moving M-point on the parameter plane can 
be used to calculate the transient response of nonlinear systems. 

As the M-point moves along the M-locus, each point defines both 
signal amplitude and all roots of the characteristic equation. 

This information can be used to determine the amplitude vs time 
relationship which is the transient response. Computations are 
based on Siljak's extension of some basic work by Krylov and 
Bogoliubov, and details are given in the following paragraphs. 

Assume that the system is second order, and that the non- 
linear element is represented by its describing function. Then 
for an initial signal amplitude A q , the transient response is de- 
fined by 

X(t) = A q e CTt cos (at + 0) (2-1) 

where a and to are both functions of the signal amplitude. 

*These assumptions restrict use of this method to systems in which 
a pair of complex roots dominate the transient response, and these 
systems must have low pass filter characteristics to justify use 
of a describing function. 



The parameter plane curves are prepared, the M-locus is super- 
imposed on them, and the describing function is used to associate 
an amplitude scale with the M-locus. Then the values of' cr(A) 
and 0>(A) may be read from the parameter plane for any X. 

The transient response of the system from any initial dis- 
placement, A q , is determined in two steps, the first of which is 
to calculate the envelope of the transient. Assuming that 0 
(in eqn. 2-1) is zero, the envelope is defined by 

X(t) = A q (2-3) 

which may be approximated over a short time interval by a straight 

line tangent to the exponential curve. Thus at t = 0, X=A q and 

from the parameter plane cr(A o )=a o is evaluated. Then X(t)=A o C CTot 
is approximated by a short straight line segment on the X vs t 

plane. This straight line is terminated at t = t^ and at t^ a 

new amplitude A^ is read from the curve. Entering the M-locus on 

the parameter plane with A^ values are obtained for and co^. 

The envelope of the transient is extended from t^ to t 2 with 

0* t 

another straight line segment defined by X = A^e 1 . This pro- 
cedure is repeated until the envelope is defined over an accept- 
able time interval. 

As a by-product of this procedure, w has been determined 


quantitatively as a function of amplitude and also as a function 
of time. Using the definition J / 



& = J W(A) dt 


(2-4) 


the phase can he determined at any t by graphical integration 
(i.e., evaluation of the area under the w(A) vs t curve). If 0 
in eqn. 2-1 is zero, then X(t) = 0 for 3? = (2n-l) (77/2). Values 
of t corresponding to $ = 90°, 270°, 450°, etc., are determined by 
graphical integration, are marked on the axis of the X vs t plane, 
and the transient response is drawn tangent to the envelope and 
intersecting the X= 0 axis at the indicated values of t. 

The above procedures are readily applied to systems with one 
nonlinearity, and correlation with simulation results is excel- 
lent. Since such applications are elementary no illustrations 
are given here, and the study is extended to systems with two 
single valued nonlinear elements. In general no other methods 
exist for predicting the transient response of systems with two 
nonlinear elements, so the results obtained here represent a 
significant advance in the state of the art. 

2.2 CLASSIFICATION OF SYSTEMS WITH TWO NONLINEARITIES 

When a system contains two nonlinear elements, N, and N 2 , 
that are single valued and are not frequency dependent, parameter 
plane representation may be used but both a and jS become functions 
of N^ and N,,. Computation of the parameter plane curves presents 
no difficulty, but determination of the M-locus may be difficult. 
As a result it is convenient to classify nonlinear systems accord- 
ing to the structural conditions which complicate the evaluation 
of the M-locus. The following classes are proposed: 


4 



CLASS 1. Identical signal excitation to both nonlinear 
elements . 

In Fig. 2 -la, the signal X is the input to both nonlinear 
elements N^ and N 2 . For every value of X corresponding values 
of N^ and N 2 are uniquely defined and are independent of fre- 
quency so evaluation of the M-locus is easy. 

CLASS 2. The input signals to the two nonlinear elements 
are related by a linear differential equation. 

In Fig. 2-lb the signal X is the input to N 2 , but the 
input to N^ is X G_^(s). Thus the input to N 2 is a function of 
amplitude only, but the input to N^ is a function of both ampli- 
tude and frequency. For a given amplitude of the signal X, the 
describing function for N 2 provides one unique value, but for 
each amplitude of X the describing function for N^ has an in- 
finite number of possible values, one for each possible value of 
frequency. As-a result the evaluation of the M-locus is consid- 
erably more difficult than for Class 1. 

CLASS 3. The input signals to the two nonlinear elements 

are related by a nonlinear differential equation. 

Fig. 2-lc illustrates this class of nonlinear systems. The 
signal X is the input to N^, but the input to N 2 is xCn^Hg^s) } 
where the brackets are intended to represent some functional re- 
lationship rather than a multiplication. Evaluation of the M- 
locus can be very difficult for such systems. 

2.3 EVALUATION OF THE M-LOCUS. THE DYNAMIC DESCRIBING FUNCTION. 

When a system with one single valued nonlinear element is 
represented on the parameter plane the M-locus is clearly a 


5 



straight, line parallel to one of the coordinate axes. Ulus the 
M-locus itself is readily found hut the amplitude scale asso- 
ciated with this locus must he evaluated.- For systems with two 
nonlinearities (especially Class 2 or 3) the path of the M-point 
on the parameter plane cannot he predicted hy inspection. It 
can he calculated, however, using the ordinary describing function 
to define the amplitude relationships. 

To justify the choice of the describing function as a tool, 
consider the fact that parameter plane predictions of limit cycles 
are defined on the basis of a single point where the M-locus in- 
terests the stability boundary. This single point defines both 
the fundamental frequency of the oscillation and also the ampli- 
tude of this fundamental component. It is clear that the loca- 
tion of the M-point represents some sort of average value of ampli- 
tude, since the instantaneous value of amplitude varies cyclically 
during a limit cycle. The describing function of a nonlinear 
element effectively averages the response of the element to a 
sinusoidal input over one cycle of operation. Thus its use is 
clearly justified when system operation is periodic and lightly 
damped. While not so clearly justified for other operating con- 
ditions it has given surprisingly accurate results and therefore 
will be used until a better technique becomes available. 

Using the describing functions of the two nonlinearities in 
a system, a family of decribing function curves are computed and 
plotted on the Oi-j S parameter plane. When these curves are super- 
imposed on the regular parameter plane curves, the M-locus can 
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be determined. The M-locus represents the curve along which the 
M-point moves when the system is in dynamic operation, and it 
consists of the locus of all points at which the describing func- 
tion curves and the parameter plane curves have common frequency 

intersections. We choose to call this curve the "Dynami c De- 

scribing Function Locus". The procedure and also a justification 
is as follows: 

a) Assume a constant amplitude, constant <0 signal at X, the 
input to one nonlinear element. Using the describing 
function compute the equivalent gain of that element? 
also compute the signal amplitude at the input to the 
second nonlinear element, and the equivalent gain of 
this second element. 

b.) The two equivalent gains evaluated in (a) determine one 
point on a describing function curve on the a-/3 plane. 
Repetition using the same value of w but different 
amplitudes at X determines a describing function curve 
for a constant co signal. 

c) Repetition of a) and b) for other values of w provides 

a family of describing function curves, each curve being 
for a designated value of oo. 


d) These curves are then superimposed on the usual* para- 


meter plane curves. The constant -co describing function 

-^Curves for constant -a and constant to are most convenient, but 
constant -£ and constant c o curves can be used if it is noted that 

“ = V 1 ^ 2 • n „ 
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curves will intersect the constant -CO parameter plane curves, 
and those intersections for which the to is the same. Define the 
Dynamic Describing Function locus. 

The nonlinear system is described by one nonlinear dif- 
ferential equation. The procedures used here effectively parti- 
tion this equation into two parts, a linear part represented by 
the parameter plane curves, and a nonlinear part represented by 
the describing function curves. Then parts are "coupled" by the 
parameters o; and /3 which are the coordinates of both plots. If 
the system is in steady state periodic motion at a given fre- 
quency the nonlinear differential equation of the system must be 
satisfied, so the linear and nonlinear partitions must be satis- 
fied at that frequency. This condition can exist only at the 
intersection of the common frequency curves. The points thus 
defined on the "Dynamic Describing Function Locus" are deter- 
mined on the basis of steady state sinusoidal operation (unforced) . 
Under transient conditions the M-point moves along some locus on 
the parameter plane, and we assume that the points on The Dynamic 
Describing Function locus apply to transient operation although 
they are determined by means of steady state sinusoidal concepts. 
Experimental results indicate that this is a good assumption. 

2.4 CALCULATED AND EXPERIMENTAL RESULTS 

In order to verify the correctness and the applicability of 
the dynamic describing function and the graphical transient re- 
sponse calculations, specific examples of each of the three gener- 
al cases of Fig. 2-1 were investigated. The details of some of 
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these examples, and the corresponding calculated results are pre- 
sented here. Simulation of the systems provided experimental re- 
sults which are also presented to permit comparison between theory 
and experiment. 

System 1 . Two nonlinear elements with identical excitation: 

The block diagram is given in Fig. 2-2. The characteristic 
equation is 

3 2 

s + 10s + ( 10N i + 10n 2 ) s + 100N 1 = 0 (2-5) 

and it is convenient to let = a, ^ = P. Fig. 2-3 gives the 
parameter plane plot (in <7— and co- curves) . Since the two non- 
linear elements have identical excitation a single dynamic de- 
scribing function curve is obtained which is independent of fre- 
quency. However, the dynamic describing function is dependent on 
the specific numerical characteristics of the nonlinearities, and 
Fig. 2-3 contains three dynamic describing function curves (dotted) 
for three different sets of characteristics in and N£. These 
three curves were chosen to illustrate different root variations. 
For curve 1 a real root becomes dominant early in the transient, 
for curves 2 and 3’ complex roots are dominant, the system being 
moderately damped for curve 2 but going to a stable limit cycle 
for curve 3. 

Calculated and analog computer results are given on Figs. 2-4, 
5,6. it is seen from Fig. 2-4 that the dominant real root condi- 
tion cannot be handled accurately with the graphical computations. 
It is not known whether the di 


repancy lies solely in the graphical 
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method which is based on complex roots, or whether the dynamic 
describing function also contributes to the errors. Research on 
this point is continuing. For the cases of Fig. 2-5 and 2-6 the 
calculated results compare well with the computer results. 

System 2 . Two nonlinear elements related to a common signal by 
a linear differential equation. 

The block diagram is given in Fig. 2-7, and the parameter 
plane curves with dynamic describing function curve shown dotted 
are given on Fig. 2-8. Fig. 2-9 gives the describing function 
grid needed to obtain the dynamic describing function curve. To 
obtain the grid of Fig.- 2-9 the point A q on Fig. 2-7 was chosen,' 
as a reference point, and a-t each value of to the amplitude of the 
(assumed) sinusoidal signal at A q was varied to obtain the vs 
N£ values for a constant to curve on Fig. 2-9. The dynamic des- 
cribing function curve on Fig. 2-8 is obtained by superimposing 
the parameter plane curves of Fig. 2-8 on the describing function 
net of Fig. 2-9 and locating intersections of constant to curves 
of the same to value. 

Limit cycle predictions of the dynamic describing function 
curve on the parameter plane agree with analog computer simulation 
results. In addition Figs. 2-10, 11,12 compare predicted transient 
response with simulation results. 

Additional checks were run using different values for the 
deadzone and saturation limits in the two nonlinearities, but 
the detailed data is not given here. In general the predicted 
and simulated results were in good agreement except when a reaX . 



root became dominant during the transient response, in which case 
the frequency of the oscillatory component was usually predicted 
with reasonable accuracy, but amplitudes were not, nor was the 
total response time due to the influence of this real root. 

The calculations and simulations were also repeated with the 
nonlinearities interchanged (i.e., in Fig. 2-7, becomes a 
saturated element and ^ a dead zone element) . Using the same 
techniques the results obtained were always in agreement with 
about the same degree of accuracy and with the shortcomings as 
previously noted. 

System 3 . Two nonlinear elements related to a common signal by 
nonlinear differential equation. 

Ihe classification described as System 3 can contain a wide 
variety of combinations of linear and nonlinear elements, of 
which the parameter plane method may be applicable to only a 
small subset. A specific system which belongs in this class is 
shown in Fig. 2-13. The characteristic equation of this system 
is 

s 3 + 3s 7 + 2s + 40KN, (N + jN, ) (2-6) 

where U 2 = N a + for the hysteretic nonlinearity, and we define 

a = |8 = P arame . ter plane equations are still appli- 

cable and the parameter plane curves can be computed. For the 
purposes of this study only C = 0 curve was calculated, and only 
the limit cycle predictions were checked. The describing function 
net is required, and in this case relates the N-^N a and pairs 

to the common signal at A on Fig. 2-13. The results of these 



computations are given on Fig. 2-14, which shows the G = 0 curve 
from the parameter plane equations and the describing function 
net for the case where K = 0.15. Only one point is defined on 
the dynamic describing function curve, and this is marked on 
the K = 0 curve at the point where the w value on the £ — 0 curve 
is the same as the value of the constant describing function curve 
passing through that point. This defines the frequency and ampli- 
tude of the limit cycle, and the results agree with simulation 
results. 

Note that a change in the value of K changes the differential 
equation of the system, thus requiring a new set of curves. Re- 
sults were obtained with other values of K and again the predic- 
tions agreed with simulation results. 

2.5 COMMENTS 

The results obtained thus far indicate that the parameter 
plane is a useful tool in predicting the stability and response 
of nonlinear systems. The accuracy available is only fair, but 
is more than adequate for many engineering applications. The 
transient response predictions - in particular for systems con- 
taining two nonlinearities, - are better than are available with 
any other method. 

/ 

The graphical presentation of the dynamic describing function 
curve on the parameter plane is potentially a valuable design 
tool. It indicates at a glance the range of variations of the 
roots, and thus permits prediction of a desired location of the 
describing function curve, which in turn implicitly defines the 

•S7 



required characteristics of the nonlinear element. Further 
research is required in this area. 

The technique becomes inaccurate when the transient re- 
sponse is influenced by more than two complex roots. Again more 
research is required to evaluate this situation. 

It is too early to assess the true value of studying non- 
linear systems on the parameter plane. Without question it does 
mahe possible many types of analyses that are not readily avail- 
able otherwise. However , the limitations of the technique are 
not clearly defined, and it obviously is important to Know under 
what conditions the. methods are not applicable, or should be 
applied with care. 
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CHAPTER III 

ASYMMETRICAL NONLINEAR OSCILLATIONS 




3.. 1 Introduction . 

In certain classes of nonlinear systems, oscillations may 
consist of a limit cycle superimposed on a constant or slow-varying 
signal. These oscillations are referred to as asymmetrical oscil - 
lations since the center of the limit cycle is shifted according 
to the corresponding value of the constant or slow-varying signal. 
In general, asymmetrical oscillations may occur when the input- 
output characteristic of the nonlinearity in the system is not 
symmetrical about the origin, or when the system is subject to 
forcing signals. When the nonlinear characteristic is asymmetric,, 
the output of the nonlinearity may contain a constant term even 
though the corresponding input is a single sinusoidal wave. If 
the nonlinear characteristic is symmetric, asymmetrical oscil- 
lations can arise whenever the system is subject to forcing input 
signals. Evidently these oscillations may take place at certain 
points of the system if both conditions are present. Before the 
analysis, of asymmetrical oscillations in the parameter plane is 
presented, the previous work and results in considering these 
oscillations and related problems are reviewed. 

It has been shown first by MacColl [3.l] that the introduc- 
tion of an external sinusoidal signal at the input to an on-off • 
servomechanism yields a system that behaves like a linear one for 
small inputs superimposed >on the ^sinusoidal 'signal. This pheno- 
mena has been later investigated under various names, such as 
"dither effect",, .7 signal 'stabilization^*, etc.' Asymmetrical non- 
linear ■ oscillations has been found by. a majority of authors as 
the most appropriate- term for the mentioned phenomena. 
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In analyzing a carrier-controlled relay servo, Lozier [3.2] 
has introduced an idea to accomplish the linearization of the re- 
lay by a. limit cycle existing in the system and without an ex- 
ternal signal. This idea has been further developed by several 
authors [3.3-3.91 and a detailed treatment of the problem has been 
given by Popov and Palitov [3.8]. On the other hand, the external 
signal application has been developed by Loeb [3.9] and Oldenburger 
with his associates [3.10-3.12], The latter introduced the name 
"signal stabilization” to indicate that the nonlinear system is 
stabilized in the- state of sustained oscillations with suffici- 
ently high frequency. The stabilization is actually a consequence 
of the linearizing effect discovered by MacColl. The concept of 
signal stabilization has been extended by Sridhar [3.13-3.14] to 
the case of a nonlinear system which has one single-valued non- 
linearity in the loop, and the stabilizing signal is a stationary 

t * • 

random process with a Gaussian distribution . and 'obeys the ergodic 
hypothesis. 

i ■ > ’ , 

The above defined problem' can be treated *by dual-input de- 
scribing functions, as proposed, by West [3.15]. This approach has 
been significantly simplified by Boyer [3.16] as outlined by 
Gibson [3.17]. A similar approach is used by Gelb and Van der Veld 
[3. 18], and significant results have been obtained by Atherton and 
others [3.19-3.20] who made a comparison of the utilized concept 
with the Tsypkin method [3.2l], 

The study of asymmetrical nonlinear oscillations has been 
extensively performed in the analysis and design of a large class 
of plant adaptive control systems. This class of system is 
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sometimes called the limit cycling adaptive systems because of 
the fact that the existing limit cycle Is used as an identifi- 
cation signal,. Some of the references on this subject are listed 
here [3.22-3.26]. A majority of the authors .proposed an external 
sinusoidal signal for identification. More recently, Gelb and 
Van der Velde [3.18] have examined to a limited extent and in 
a quantitative manner the properties of self -oscillating adaptive 
systems which have several advantages over the external adapta- 
tion, such as simplicity, cost, reliability, etc. The following 
analysis of asymmetrical nonlinear oscillations in the parameter 
plane can be applied directly to self-oscillating adaptive sys- 
tems. 

In the following developments, the asymmetrical nonlinear 
oscillations are analyzed in the parameter plane [3.27], The 
control systems with asymmetrical nonlinear characteristics are 
considered to determine stability and sustained oscillations. 

The same type of oscillations is investigated in nonlinear con- 
trol systems subject to constant reference and perturbing input 
signals. The procedure is further extended to the analysis of 
systems with slow-varying input signals.' In this case, it is 
shown .how a nonlinear characteristic can be. modified for the 
slow-varying signals. , The presented analysis is performed with 
respect to both input signals and the values of adjustable sys- 
tem parameters. The analysis procedure is illustrated by examples 
in which multiloop feedback structures with several adjustable 
parameters are considered. In addition, various nonlinear 
characteristics are used in either the forward or the feedback 
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path. The obtained results are checked by computer simulations 
which indicate a sufficient accuracy of the presented procedure. 
3.2 Basic Developments 

Consider a nonlinear system described by the nonlinear dif- 
ferential equation 

B ( s.) x + C(s) F(x,sx) = H(s)f, s s ^ (3,. 1) 

where B(s), C(s), and H(s.) are polynomials in s and the degree of 
the polynomial B(s) is -greater than the degree of the polynomials. 
C(s) and H(s) . The function F(x,sx) describes the nonlinearity. 
Function f = f(t) is a forcing signal, which may be 'either a 
reference input or a perturbing signal, and it is assumed to be 
a constant or a slowly-varying function of time. 

As a first approximation, the steady-state, solution x = x(t) 
of equation 3 . 1 which' 'represents the input to the nonlinearity, 
is assumed to be 

x = x° + x* (3.2) 

where x° = x° (t) is either a slowly-varying function of time or 
is constant, and x*, which is 

x* = A sin 0, 0‘ = Qt + 0, (3.3) 


represents the periodic component -of the solution x(t). Since 
0' in (3.3) :merely corresponds to a shift in t, one can put 0=0 
and use x = A sin ^t. 

The forcing function f(t) is considered as a slowly- varying 
function of time if it can be assumed approximately as constant 
over any cycle of the periodic component x*; i.e., 

|f (t+T) - f (t) | < < if(t) | (3.4) 
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where the period T = 277/£2. In the frequency domain, equation 3.4 
means that the frequency £2 of the periodic component x is much 
greater (practically ten times or more) than the highest frequency 
of the slowly- varying component x°. In this case, no harmonic 
relation between the components x and x nonlinear system sub- 
ject to forcing signals, such as jump-resonance,' generation of 
subharmonics, etc., cannot take place. The forced nonlinear os- 
cillations for which the condition (3.4). is not satisfied neces- 
sarily, are considered in other works. 

Under the condition (3.4), the values of x°. A, and £2,. which 
appear in' the solution x = x° + A sin £2t, are slowly-varying 
.quantities in time. This enables the extension of the conven- 
tional harmonic linearization in which the describing function 
is defined for the signal x — x° + x as an input to the non- 
linear element. ‘Thus, ,the nonlinear function ^F(x>sx) is- approxi- 

' ’ ‘ ' ' * ‘ * / ’ - y ** 

mately expressed by the first terms of the Fourier series as 

i * 

F(x, sx) •= F° + N^x + "jjjysx (3.5) 

. 1 * ' ' 

where 

1 r 2lT 

F° ~ j F(x° + A sin 0,, AO- cos 0)'d r 0 (3.6a) 

0 o 

l P 21T 

Nf = ^ j F(x° + A sin 0, A £2- cos 0.) sin 0 d0 (3.6b) 
o • 

1 r 27r 

N 2 = j F(x° + A sin 0, A £2 cos 0)cos".0 -d0 (‘3.6c) 

' o 

and 0 = £2 t . 

As can be seen from- equations 3.-5 and 3.6a, the component 
F° of the output of the nonlinearity F(x, sx) is not considered: 
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zero as was the case in the analysis of symmetrical nonlinear 
oscillations presented in the previous chapter. This results 
from the fact that either the nonlinear function F(x, sx) is not 
symmetric or the system is subject to an external input signal,, 
or that both facts are present in the system. 

According to equations 3.6, all coefficients E°, N-^ and N 2 
.are generally functions of x°. A, and i.e., 

F°+F°(x°,A,Q) , N 1 - N-^x^A,^, N 2 =N 2 (x°,A,0), (3.7) 

For a majority of the nonlinear .functions F(x>sx), encountered in 
practical applications, the above functions (,3-. 7) are obtained 
once and for all,. 

By applying the linearization of the function F(x,sx) given 
in equation 3.5, the solution x = x° + x of (3.1) can be ob- 
tained by considering the following linearized differential 
equation 

B(s) (x°+x*) + C(s) (F^N-jX* -+ sx*) = H(s)f (3.8) 

instead of equation 3.1. If x°, A, and are slowly-varying 
functions' of time as a consequence of the same property associated 
with the forcing function f, equation 3.8 can be rewritten as two 
simultaneous equations corresponding to the slowly- varying sig- 

Q ^ 

nal x- and the periodic .signal x as follows: 

B(s,)x° + C(s)F° = H:(s,)f ^ (3.9a) 

Btsjx + e(s) (•N 1 x* + SX*). = 0 (,3.9b). 

' -i 1 f \ . 

Equations 3.9, however, cannot be solved independently since they 
are related to each other by ‘the nonlinear •equations’ ’3'.?. This 
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fact indicates that the applied linearization preserves the 
essential feature of nonlinear systems and that the superposi- 
tion principle from linear analysis is not valid. 

An analytical solution of equations 3-. 9 is difficult to 
obtain since F° in- (.3 .9a) is usually a trancendental function 
with respect to x°. A graphical procedure is presented for 
solving equations 3.9 in the parameter plane. A necessary con- 
dition for equation 3.1 to have a solution x(t) close to 3. '2 is * 
that the characteristic equation 

Bts) + C(s) (N x t —;) = 0 (3.. 10) 

corresponding to the linearized differential equation 3.9b,, have 
a pure imaginary root s = jO. 

By using the parameter plane approach, equation 3.10 can be 
solved for a and as 

a-=! Oi'(Q) (3.11) 

0 = 3 (^ 

where a and $ are N, and N~ or some other system adjustable 

. • - i “ 1 

* * ■ i i , . « ‘ 

parameter. Equations 3.11 represent the S = 0 (or £ = 0) curve 
for which s « jO. The £ = 0 curve determines the stable region 
in the- plane in the usual manner,. After the stable region is 
found, the loci of points M(a, jS) are plotted according to- the 
variations of & and/or £ representing N^, and/or N 2 - The M loci 
incorporates the additional variable x°, and a family of the 
loci should be constructed for different values of x°. Then the 
stability of the nonlinear -system is determined by the relative 
location of the £ curve and the M- loci and the limit cycles are 
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found at their intersections. The stability of the limit cycles 
is determined in the usual manner. This part of the solution 
process will be best described by the examples that follow. 

The presence of a limit -cycle in the system can modify the 
nonlinear characteristic for the slowly -varying input signal. 

In order to determine the modified characteristic, the intersec- 
tions of the £ = 0 curve and the M loci are considered to evalu- 
ate the amplitude A and, the frequency £2 of the limit cycle as. 
functions of the slowly-varying component x°; i.e., 

A = A (x°) , 0 = 0(x°) (3.12) 

These functions, when -substituted into the function F°(x°,A,£2) 
yield the modified nonlinear characteristic for the slowly- 
varying signal 

F° = *M'x°) <(3.13) 

The function ib (x°) is continuous 1 , in a limited range of x°“, ' which 
indicates the smoothing effect due to the presence of the limit 
cycle. 

Substitution of equation 3.13 into equation 3., 8a gives 

* 1 ' . * 

B(s)x° + C (s) *(x°) = H ( s )• f (3.14) 

Equation 3.. 14 is a nonlinear differential equation in x°, which 
can be solved graphically for x° after the function $ (x°) is 
obtained. This, in turn, yields ‘ the related values of the 
functions A'(x ) and £2(x°) of equations 3.12, and the solution 
x = x° + A sin Ofc is thereby determined. 

The function ib (x°) is a continuous function of x° and it can 
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be assumed approximately linear for small variations of x°. Then 
the stability problem related to equation 3.14 can be solved by 
known linear methods. If it is regarded as a nonlinear function 
of x 0 , it can be linearized by harmonic linearization and the 
results of the previous chapter can be applied. 

It should be noted here that the same parameter plane pro- 
cedure can be used when the right side of equation 3.1 has more 

than one forcing function; i.e., the right-hand .side is expressed 
r 

by £ H. (s.)f, . The solution x, however, must be found by con- 
i=l 1 1 

sidering all existing inputs simultaneously since the super- 
position principle of linear analysis is not valid. Further- 
more, if the polynomial H(s) of equation 3-1 can be factored in 
the form sH^(s), the procedure applied to the case in which the 
rate sf of the function f is .considered as a slowly- varying sig- 
nal.; i.e.,, |sf(t+T) - sf (t) |. 

The presented graphical procedure can be extended to non- 
linear control systems with two nonlinear functions F '(s) and 
■F 2 (x) , whereby the following nonlinear differential equation is 
investigated : 

B { s) x + C(s) F-^x) + D(s) F 2 (x) = H(s) f . (3.15) 

In this case, a procedure similar to that given in Section 

o * 

can be extended to determine the solution x = x + x . 

The general procedure outlined in this section is modified 
depending on the actual problem involved. These t problems ' may be 
divided into three major groups: asymmetrical nonlinearities; 
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constant forcing signals; and slow- varying signals. In the 
following, each group is considered separately. 

3,. 3 Asymmetrical Nonlinearities . 

In an autonomous’ nonlinear system, which is described by the 
differential equation 3.1 and where f = 0, the asymmetrical os- 
cillations may occur whenever the function F(x,sx) is not sym- 
metrical to the origin. Then, under the conditions discussed ’in 
the previous section, the system may be described by equations 
3.9 which has the form 

B(pj x° + C(o)F° = O’ (3.16a) 

CB(s) + C (s) + -q-s-) Jx = 0 (3.16b) 


In equation 3. ,16a-, which corresponds to, equation 3.9a, there is 
no forcing slowly-varying function (f = 0) , and in the steady- 
state solution x = x + x , the x p is constant and hence s is 
replaced by zero in B(s) and C-(s) . 

In practical situations, B(o) or C(p). can be zero. Also, 
the nonlinearity in the system is often described’ by a single- 
valued function F(x) and' N 2 =Q. Thus, an adjustable parameter 
appearing in B(s) or C(s) can be chosen as one of the axes in 
the parameter ctfi plane, while the other axes is related to the 
describing function coefficient N-^. Some of these situations are 
discussed in the following examples. 

Consider a feedback control system with the block diagram 
of Fig. 3.1 in which the transfer functions are 

Gl (s)=Ki, G 2 (s)= r|_, g 3 = G^CsHK^s. .(3.17) 
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The nonlinearity n has the form shown in the upper left corner 
of Fig-. 3.2'. 

, * • ■ ! 

Equations 3.16, for the system under investigation, have 
the form 

F° =■ 0 (3.18a) 

Cs (s+1) (s+2) + [■K 2 K_ 1 s (s+2) + K^K^N^x =0 (3.18b) 
where, according to the function F(x) of Fig. 3. .2 and equations 


3.6,, one has 



(3.19a) 

(3.19b) 


N 2 s o (3 . 19c) 

and x = x(t') is the input signal to the nonlinearity n as indicated 
in Fig. 3.1. 

The characteristic equation of equation 3.18b is 

s (s+1) (.s+2) + [K 2 ,K_ 1 s(s+2.) + = 0 (3.20'). ' 

By denoting K 2 K_^N^ = & and = j3, the £ = 0 curve is 

obtained as 

a = -|(^ » 2) (3.21) 

/S = + 4) 

and the stable region is determined in the plane in the usual 
fashion as shown in Fig. 3.2. 

From equations 3. 18a and 3.19a, one obtains 

x° = ,A cos — ' (3.22) 

1+m 




and of equation 6.19b becomes 


N 


1 



sin 


TT 

1+m 
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(3.23) 


By using equation 3.23 and the expressions a — K 2 K_^liL, 

$ = three M loci (a) , (b)., and (c) , are drawn in Fig. 

3.2. They correspond to the parameter values m = 0.5, c = 1, 

K 2 = 1 and (a) K-jK^ K_ ± •= 0.125; (b) = 8.39, K_ x = 0.28? 

(c)- K-^K^ — 26, K_^ = 1.75. The stable asymmetrical oscillations 
are found at the point and M 2 where the M loci (a) and (b) 
intersect the £ = 0 curve. The amplitudes' of the oscillations are 
approximately A^ = 0.85 and A 2 = 0.8, which is read from the 
M loci (a) and (b) at the intersections and M 2 . The corres- 
ponding frequencies = 1..5 and ^ 2 = 1.6 are indicated on, the 
f - 0 curve. The related values of x° in the solution x^x+ffil’^in 

, i 1 1 » , 

is calculated for each point: and M 2 using equation 3.22, namely, 

x° = -0.42 and x° = -0„39. ) 

In Fig. 3.3, the solution x, = 0.42 + ,0.85 sin 1.5t for the 

s’ - ■ ’ • • 

case (a) is shown as, obtained by a digital' computer simulation. 

The calculated results are sufficiently close to that obtained 
by the simulation. From Fig. 3.3, it can be seen that .an initial 
condition x^(o) = 4.25 is used and the variable x^(t) approached 
a stable limit cycle. That the limit cycle is stable and will be 
reached by x^ (t) starting from x^(0) = 4.25 can be > concluded from 
the relative location of the £ = 0 curve and the M locus (a) , as 
explained in- the preceding chapter on the symmetrical oscilla- 
tions. . The additional component x° of the solution x(t) does 
not ali£er the stability analysis of the oscillations. 
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An analog computer .simulation of the case (b) gives the 
solution x 2 ■ =' -0’.39 + 0..8,-sin’ list as shown in Fig. 3.4. A suf- 
ficient accuracy ,is indicated. The initial condition x„(o) = 0 

• ‘ ; * . * ‘ . * . ' ' * h Z 

and x 2 (t) reached a limit cycle. This could he concluded from 
Fig. 3.2 as previously noted. 

It is of particular interest to consider the case (c) of 
Fig. 3.2.. The M locus (c) is tangent to the = 0 curve and cor- 
responds to the ratio a/| 3 = = 14.8. If this ratio is 

higher than 14.8, then there is a limit cycle as shown by cases 
(a) and (b) . On the other handy if this ratio is less than 14.8, 
the entire M locus is situated in the stable region and the cor- 
responding system is always stable. The tangent case (c) : 
m = 0.5, c = 1, K 2 = 1, =26, K_^ =1.75, is simulated on • 

a digital computer and the obtained .solution x^ (t) is shown in 
Fig. 3 .*5,1 which indicates that the system, is stable. 

3.4 Constant Forcing Signals 

When the forcing signal at certain points of a nonlinear 
system is constant, the solution x = x° + A sin £2t(if it exists) 
will have x°. A, and Q as constant values. To determine these 
values, note that the equations to solve in the presence of a 
constant forcing signal f° have the form 

B,(o) x° + C ( O') F° = H(o) f° (3.24a) 


■N ? * 

[b(s) 4- C(s)N^ + -Q- s) ]x = 0 (3.24b), 

In general B(o), C(o)., and H(o) are constants different 
from zero, and the solution procedure is somewhat more compli- 
cated to perform than in the previous section where the right 
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side of equation 3.24a was zero.. 

To illustrate the solution procedure, consider a nonlinear 
feedback system with the block diagram of Fig. 3.6 and the trans- 
fer functions 


G h s > 


0.2s +0.8S+1 


t 


0. 5 (s+1) K -1 

0.2s+l ' -1 T_ 1 s+1 


(3.25) 

The nonlinearity n is given in Fig. 3.7a. The input to the system 
is a perturbation signal f = f(t) which is -related to the signal 
x = x ( t.) and c = c(t) of Fig. 3.6 as 

(0.2s+l)c = 0.5 (s+1) x-f (3.26) 

If the perturbation signal is f(t) = f° = const., equations 3.24 
have the form 

x° + K_ 1 F° = f° (3.27a) 

(0.04s 4 + 0.36s 3 + 2s 2, + 2s.) T _’+() ..4s+2)K _ 1 N 1 + 


" + 0.04s 3 V 0.36s 2 ' + 2s + 2 = 0 (3.27b) 

where equation 3.27b represents' 'the characteristic equation of the 

linearized equation 3.24b. By substituting T_^ = &. and K 

the parameter plane diagram is plotted in Fig. 3-7b according to 

the parameter plane equations 

a = 0..64S? + 3.2 

. 0.016S# - O.Ogsf -4 


fi = 


0.016O 6 - 0.03d 1 + 2.56 C? 
0.016S# - 0.08ft 2 - 4 



+ 4 


• (3.28) 
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The variation of the M point due to the function - N^(x°, A) 
given as 

,, * k , . D-x° , . D+x° 

N-^ = k - —(arc sxn — — f- arc sin ~ — + 

+ ^ A - (^) ? + ^ A - (^> 2 , A*D + |x°| 

(3.29) 

is plotted in Fig. 3.7c. (The expression (3.29)', corresponds to 

the nonlinearity of Fig. 3.7a) . In order to find a solution 
o * 

x = x + x of equations 3.27., the parameter k is assumed equal 
to one, and the function F°(x°, A) is plotted in Fig. 3.7d by 
using 

F o _ 

. kp^ • D-x° . D+x° v 

+ — LD(arc sm — — — - arc sin — ) - 


AT^ 


, D-x° ' 2 
' A ' 


-A 


<^i 2 + ta° + 


o,^ . D-x° , D+x° a i . | o i 

x (arc sm “ — + arc sin- — ) J, A^D + |x | 


A 


(3.30) 


For T_p = 0 ..04', ‘ the' point Mp«'(i);..04 j 14.3) corresponds to a 

, o * _ 

solution x = x + x which will have U =• 12 rad/sec as indicated 

, \ t * . ' , * * 

on the curve C = 0‘." If K_^ =,20,"' from M,^ it follows that 
Np = = 0.715. This value of determines the relation- 

ship between the values of x° and A for a possible solution x. 
This relationship, expressed as a function A = A(x°), can be 
graphically obtained from the diagram = N^(x°, A) by plotting 
the straight line PpP 2 corresponding to the value = 0.715. 

The function A = A(x°) represents the solution of equation 
3.27b only. The pair of values (x°. A), .which tenter into the \- 
actual solution of equation 3.27,. is replotted on the diagram 












•} t. 




£ig._3 .5— r-JNonl j near-characteristic,: 



3-25 


F P = F°'(x 0 , A) of Fig.. 3.7-d into the curve P^P 2 . Suppose that 
the constant perturbing signal has a value of f° = 11.75;, then 
equation 3.27a determines the straight line f°' = 11.75 plotted 
in the diagram F° = F°'(x°, A) of Fig. 3V7d. The intersection R 
of that straight line and the curve P^P 2 9i ves the pair (x°, A) 
of the solution x(t) which satisfies -equation 3.27 simultaneously. 
At this point R, the values are x°/D = 1.35 and A/D =1. The 
same values are obtained at the point Q on the diagram N-^ = N 1 (x° / A) 
and the solution x = x° +' A sin- fit' of equations 3.27 is found. 

If D = 1, it is x = 1.3)5 + sin* 12t: 1 Note ,thht the same solution 

» 'i . ; . , , , . 

* * * * * V * * 

is obtained if the point M 9 of Fig. 3.7.b is considered save that 
the frequency O is Tower (approximately fi =; 5.5 rad/sec,). 

Simpler situations may occur if one- of the values B(o) or 
C(o) is zero. To -illustrate, consider the nonlinear system of 
Fig. 3.8. ' The transfer functions are 


G-j.Cs) 


K x , G 2 (s) 


K 2 

s(s+l)' G 3 (s) ' 



G_i(s) 



(3 . 31) 


a'nd the nonlinearity n in .the system is given by the function 
F(x) of Fig. 3.9. The input to the system is the reference con- 
stant input signal r(t) = r°. 


The nonlinear differential equation describing the above 
system is 

[s (s+1) (s+2)+K 1 K 2 K 3 ]x+K 2 K_ 1 s(s+2)F(x) = -^^(3+2)^ 

(3.32) 


which may be rewritten according to equations 3. .24 as 
K i K 2 K 3 x ° = 2r ° 



(3.33a) 
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Cs (s-fl) (s+2) + K 1 K 2 K 3 + K 2 K_ 1 s (s+2)N 1 ]x = 0 (3.33b) 

The characteristic equation of the equation 3.33b is evidently 


S'('Sfl)' (s+2) + K 1 K 2 K 3 + K 2 K_ 1 s (s+2)N 1 = 0. 
By denoting 

* “ K 3 


(3.34) 


(3.. 35) 


the parameter plane diagram is plotted in Fig. 3.10 in the usual 
fashion. The function N^ = N^(A,x°), which appears as a varia- 
tion of 0L in the point M(a; j3) is plotted in Fig. 3. 11 by using 
general formula 3.6b. 

From equation 3.33a, one can derive the following relationship 
between the input r°, the constant term x°, and the parameter 
0 = k , 

2r° 

S|8 = (3.36) 

x°/S 

where S is the parameter of the non-linearity ; F (x) , of Fig. 3.9,. The 

i ’ * i * * 

function S$ given in (3.36) is plotted in Fig. 3.10. 

f T 

' < , * ' ’ V t » - 

Now, by using ,Fig. i 3.10 and 3.11, it is possible to determine 

• , t • 

the sustained oscillations and their stability for various values 
of system parameters K^, K 2 , K 3 , K_-^, S, k, and the input r°. For 
example, if = 1, K 2 = 10, K 3 = 1.75, K_^ =1/ S = 1, k = 1, and 
r° = 1.1, then the solution of equation 3.33 is determined by the 
values x° = 1.2, A - 0.3, and O = 2.1 rad/sec to be approximately 

x = 1.2 + 0.3 sin 2_.lt (3.37) 

For a given value of ,6 = K 3 = 1.75, r° = 1.1, and S = 1, the value 
of x° = 1.2 is read from the left part of Fig. 3.10. Then the 
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value of K^K 2 0 = 17.5 determines the point M(1.2? 17.5) on the 
C = 0 curve where SI =*2.1 rad/sec.; 1 , 'At this point, K1K .a = 1.2 

which gives N. = a = 0.12. Fig. 3. 11, is used ,to evaluate the 

* . - 0 , ; • ‘ a ' I . . * 

amplitude A, = 0..3 from the curve x /S = 1.2. The value A ~ 0.3 

is read directly from the diagram ‘N T (A, 1 x°) , of Fig.l 3.. 11, since 

: ■ • v > , • - > ’ ; „ 

K 

K = S = 1 are the parameters of the given nonlinearity in Fig.. 

3.9. 


The solution (3.37) is stable since an increase in the 
amplitude A causes the point M to move into the stable ‘ region,; 
while a decrease in the amplitude A places the' point M inside 
the unstable region of the parameter plane ('Fig. 3.10). It is 
of interest' to note that if the - produce K.K 2 ^ where $ = is 
such that it is less than 6.4, the system is always stable since 
there is no intersections of the variation of the. M, point and the 
C = 0 curve. 

The above solution (3.37) is checked’ by computer simulation 
to obtain the curve on Fig. 3.12. The accuracy of the calculated 
solution is sufficiently high and, for calculated values of x°. A, 
and Si, is approximately 10%. On the other hand, the computer 
solution indicates a distortion of the assumed solution 
x = x G + A sin (i which is due to the higher harmonics present 
in the actual solution. 

3 . 5 Slowly-varying Signals 

In this section, the problem, of linearizing a nonlinear 
system by a high-frequency limit cycle is considered in more 
detail. The objective is to determine the conditions under which 
















Parameter. .Diane _diacrain__ 
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such a linearization is possible and then to construct the 
linearized characteristic of the nonlinearity. This lineari- 
zation has several practical aspects discussed in Section 3.1, 
which are based upon a general property of the linearized system 
that, for a limited magnitude of the reference signal, behaves 
like a linear system. Therefore, results of the nonlinearities, 
such as dead-zone, hysteresis, backlash, etc.., are eliminated. 

The procedure to achieve this will be best illustrated in the 
following examples. 

Consider the system on Pig. 3.13 with the transfer functions 


K, 


K, 


(s) — K, G 2 (s) — — ^ 


s +0. 8s+l 


' G 3 (s) s (s+1) 


, G_ x (s) 1 = K_ 1 


(3 .38) 


and the nonlinearity n as shown in Pig.. 3.14. The input to the 
system is a slowly-varying reference signal r = r (t) . 

The equation which describes the system is 

[s ( s+1) (s 2 +0.8s+l)+K 2 K_ 1 s (s+l) ]x + KjK^KjFfx) = K^s (s+1) r 

(3.39) 

where the signal x = x(t)' is the 'input , to. the’ ponl inear ity. Equa- 
tion 3.39 can be rewritten in terms of equations 3.9 as 

Cs(s+1) (s +0.8s+l)+K^K^s (s+l.)]x°+K 1 K 2 K 3 F° = K^s (s+1) r 

[ s (s+1 ) ( s f0.8s+l)+K 2 K_ 1 s (s+1) l+K^^^ x = 0 (3.40) 


The characteristic equation of the second equation 3.40 is 

s(s+l) (s 2 +0.8s+l)+K 2 K_ 1 s(s+l) + K 1 K 2 K 3 N 1 = 0 (3.41) 

Substituting K 2 K_-^ = esc, K^K 2 K 3 N^ = /3, and s = j‘!T2 into equation 
3.41, one obtains the parameter plane equations of the £ = 0 curve 
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as 

a = 1.8 & - 1 

(3.42) 

j3 = 0.8 0{O + 1) . 

The curve C = 0 is plotted in Fig. 3.14. • The variations 
of the M point are plotted also in Fig. 3.14 according to 

N 1 = If ^ ~ { | )2 ' A * U° 1 (3.43) 

The system parameters = 1, K 2 = 12.5, K 3 = 10, K_^ = 1 
result in the point M-j, (12 . 5 ; 45). If c = 1, this point gives 
= ^/K^K 2 K 3 = 0.36, and the straight line P' 1 P 2 pl’ otte< 3' on 
the diagram of function (x°, A) . After the diagram 

F° = F°(x°, A) is plotted in Fig. 3.15 using 

F° = -^ arc sin 4^ , A a l x °l (3.44) 


the replotting of the straight line P-^P 2 on diagram F°(x°, A) 

yields the function (x°) of Fig. 3.15. The replotting procedure 
is the same as that used in the ‘ previous section; i.e., for each 
pair of values (.x°,A) read on the straight line P^P 2 , the cor- 
responding pair exists in the diagram F°(x , A ) „ which determines 

one point on the c'urve $(x°)V 

' « * 

Function ii(x°) of Fig. 3.15 is smooth and represents the non- 
linearity for the slowly- varying signal x Q . The shape of ip (x°) 
explains the smoothing effect of the high frequency limit cycle 
which has a slowly-varying amplitude, the value of which is 
located between the points and Q 2 on the A axis of Fig. 3.14. 
The frequency ^ is approximately constant and has the value 
O -- 2.7 rad/sec. According to ip(x°), the smoothing effect of the 
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limit cycle is present under the' condition that |x w |^ 2',..25-. For 
small values of x°, it is possible to consider !/) (x°) = Kx° where 
K = const. Then the stability of the system with respect to 
slowly-varying signals may be investigated by well-known- linear 
methods outlined in. Chapter II. t In the specific example, the 

l 

equation of interest is 

s ( s+1) (s +0 . 8s+l) ' -i- K„K - s (s.+l) 1 + ; k' K-K ‘K. = 0 (:3.45) 

Z 1 Z J 

Finally, it is to be noted that for the smoothing effect to take 
place, the amplitude A should be A ^ |x°|, as stated in equations 
3.43 and 3.44. 

The results of the above analysis are checked by simulating 
the system on an analog computer. Three cases are considered. 

In Fig. 3.16, the input to the nonlinearity x = x° + A sin Ot 
and the system output x = x(t) are shown when the input signal 
is r-= sin O.lt. The obtained computer solution agrees with the 
prediction. The output c(t) exhibits a smaller amplitude limit 
cycle with the same frequency. When the input amplitude is in- 
creased five times, the diagram of Fig. 3.17 Is obtained. This 
change increased. x°, but the amplitude A remained almost the 
same,. The frequency Cl did not change. Similar results occurred 
when the input amplitude increased ten times except that the 
amplitude A became slightly smaller, which agrees with the dia- 
gram of Fig. 3.14. The third case is given in Fig. 3.18. It 
should be noted from these’ computer ‘solutions that the output 
signal c (t) represents the input signal r(t) except for the 
superimposed limit cycle. It can be eliminated by introducing 
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sufficient filtering in the block (-s) of the system of Fig. 3.13, 
or by readjusting the system parameters to obtain a higher fre- 
quency limit cycle. 

If the values of the system parameters are chosen so that the 
operating point is M 2 (21.2? 120,) • of Fig. 3.14, the frequency of 
the limit cycle becomes higher. However, the corresponding range 
of variations of x° is decreased to |x°| <0.7, together with the 
range of the amplitude A which is between Q 3 and Q^. This indi- 
cates that the presented procedure' is convenient to apply when 
the system parameters and operating conditions are changed. 

, . . * * # 1 t \ » * " t i ,1 

If the nonlinearity h is changed in the system of Fig. 3.13 
by introducing a considerable dead Uorie, D/ a' diagram of Fig. 3.19 
is obtained. The variation of the M point is calculated by using 
equation 3.6b for the given nonlinearity of Fig. 3.19. Two cases 
should be considered separately; i.e., 

N 1 = + A &|x°j+D (3.46a) 

N 1 " If Vl-C 2 ^) 2 , [x°|-D < A <’jx°f.+D (3.46b) 

and the diagram N-^x 0 , A) is shown in Fig. 3.19. By using equa- 
tion 3.6a, the corresponding diagram, F.(x°, A) of Fig. 3.20 is 
plotted according to 

F° = § (arc sin — A ~ - + arc sin , A s=|x°| + D (3.. 47a) 


= £ (£- + ' 
7 T ^2. •' 


arc sin 


x -D 


A 


) sin 


o 

x , 


■xr 


-D- £<A < x +D 


(3.47b) 
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If the points and Mj are chosen in Fig, 3.19 as operating 
points, the replotting of the straight lines and PjP* results 

in the two linearized characteristics a and b of Fig. 3.20, re- 
spectively. They are constructed for the values of nonlinear 
parameters c = D = 1. As can be seen from Fig. 3.20 the dead 
zone is eliminated as far as the slowly— varying signals are 
concerned. For this to take place, it is necessary to choose 
operating conditions such that equation 3.47a is valid. This 
means that the amplitude A of the limit cycle must be greater 
than !x° l+D. Otherwise the linearized characteristic ii (x°) does 
not go to zero when x° = 0 since F° does not go to zero for 
x° = 0. This is indicated in Fig. 3.20 whereby F° = 0 for x° = 0 
and the dead zone is eliminated. 

By the outlined technique, it is possible to eliminate the 
hysteresis and backlash in systems with multi-valued nonlineari- 
ties. The linearization yields a single-valued function t (x°) 
which is linear in a certain limited range of values of the vari- 
able x° about the origin. To illustrate this, consider a non- 
linear system with the block diagram of Fig. 3.21 and the trans- 
fer functions 


The nonlinear function F(x) of the nonlinearity n is given in 
Fig. 3.22. 

The equation describing the system is 



(3.48) 


s (s+1) (s+2) s + (K,s ,)K.F(x) = 0 


'-1 - l '*! 


(3.49) 



After harmonic linearization of 3.49, the corresponding 


characteristic equation is 
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If = 50, K_ 1 


a = N, 


s (s+1) (s+2) + K-^K^s+l) (N + -j^s) = 0 
1 


(3.50) 


(3.51) 


and s = j 0, one obtains the C = 0 curve as 




S = -i-fi 
p 2!> * 


(3.52) 

The curve is plotted in Fig. 3.22. On the same plot, the varia- 
tion of the point M(N^; N 2 ) is constructed according to 

(v 


Nj = — (Jl- 


D - x ^2 + / j 7 , 0 + x ";2 


-)‘ 


) ) A = D + 


N 2 = 


4cD 

»TA 2 


(3.53) 


and the nonlinearity F(x) of Fig. 3.22 for which c = D = 1. From 
the intersections of the C = 0 curve and the variation of the M 
point, one can determine the amplitude A and the frequency 0 
as function of x°; i.e., 

A = A(x°) 

0 = (x°) (3.54) 


Then, by using the expression 
F° - f (^c sin ^ 


arc sin ~ ) , A * D + |x°| 


(3.55) 


for c = D = 1, a family of curves with constant amplitude A is 
plotted on Fig. 3.23. If the first equation 3.54 is mapped onto 


/03 






































Fie. ^-25 - Comnuter solution for x a 2„6 sin 4,8t 
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the family of constant amplitude A, the function- 0 (x°) is 
obtained as shown in Fig. 3.23. The function $ (x°) as a single- 
valued function of x°', which is linear in the range- o^jx°|^ 2.4. 

For an input r(t) - 5 sin 0.5t, the computer solution is 
shown in Fig. 3.24. The amplitude A and the frequency O of the 
limit cycle are slowly-varying quantities according to equations 
3.54 and the slowly- varying variable x°. Their average values, 
however, are close to that which can be predicted from the parameter 
plane diagram of Fig. 3.22;- i.e., A = 2.8 and 0~'4.5 rad/sec. 

This can be concluded from the diagram (a) of Fig. 3.24. On- 
the diagram fb.) , the output signal c (t) is shown whereby the 
limit cycle is largely attenuated by the block G^(s-) of Fig. 3.21. 
The low-frequency component in the signal c(t) represents the 
input r(t), = 5 sin 0.5t at the output of the system. 

Of course, if the input r(t) is not present, the system- 
will exhibit a limit cycle which can be determined from the inter- 
section of the M locus x° = 0 and the £ = 0 curve on -Fig. 3.22 
as x = A sin t, A = 2.6-, O = 4.8. This is checked by the analog 
computer simulation and the obtained solution is shown on Fig,. 3.25. 
3 . 6 Conclusion 

The parameter plane method has been used to indicate ex- 
istence of asymmetrical oscillations' in nonlinear control systems. 

A procedure has been developed determine j the '’oscillations for 
different values of system parameters and input . signals. It has 

* , » , » . I 

been shown how a limit cycle -can modify the nonlinear character- 
istic for slowly-varying signals. This modification may be of 
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importance when a high-accuracy control system has to be de- 
signed in. the presence of nonlinearities with excessive dead 
zone, hysteresis, backlash, etc. The design- technique can be 
directly applied to a large class of plant-adaptive control sys- 
tems where a sinusoidal signal is used as an identification sig- 
nal. 

In a future study, the technique may be extended to the in- 
vestigation of transient asymmetrical oscillations. Thus, to 
study how these oscillations are established after certain ampli- 
tude perturbation, this study should be largely based upon the 
material presented in the following chapter. , - 

It may also be shown [16, 17 1 that the presented analysis can 
be extended to the case when the signal superimposed on a sinusoid 
is not only a constant or slowly- varying sinusoid, but also when 
the additional .signal is described as. a Gaussian process, pro- 
vided that the amplitude or standard deviation of the additional 
signal is of no consequence in the analysis. This further gener- 
ates the idea of applying the dual-input describing function 
[15,17] along with, the parameter plane method, and investigates 
the case when the input to a nonlinearity of the system is a 
combination of two similar sinusoidal signals. 
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Classical techniques for analysis and design of dynamic 
systems are largely restricted to cases m which only one para- _ 
meter of the system is adjustable. As a consequence complex 
systems cannot be treated adequately with classical techniques. 
Algebraic methods, as developed m NASA CR-617 , are capable of 
treating systems m which two parameters are adjustable, and 
thus permit analysis and synthesis of systems which are too 
complex for treatment with classical methods. 

The treatment of algebraic methods presented in CR-617 
develops the fundamental theoretical basis for the coefficient 
plane and parameter plane methods. It also applies these 
methods to basic problems such as stability analysis, cascade 
compensation of..systems, and related - topics. The applications 
indicated in CR-617 are rather elementary, i.e., the problems 
considered illustrated the procedures to be used but were not 
very complex problems. This report is based on the findings of 
CR-617, and extends the applications of the algebraic methods to 
problems of a more complex nature. 

When cascade compensation is used in a feedback control 
system, more than one filter section may be required to achieve 
desired performance. Frequency response methods involving trial 
and error are often used, but parameter plane methods permit 
analysis and design without trial and error if it is permissible 

^Algebraic Methods for Dynamic Systems by G. J. Thaler, D. D. 
Siljak and R. C. Dorf, Nasa Contractor Report NASA CR-617, 

Nov. , 1966. 



to use two identical filter sections. This problem is treated 
in Chapter I of this report. The applicable parameter plane 
equations are derived and a digital computer program based on 
these equations is presented. The program is used to study the 
effects of compensation on several systems. 

Chapters II and III are concerned with nonlinear systems. 
Conventional methods such as frequency domain analysis of sys- 
tems with the Describing function have proven useful when the 
system contains only one nonlinearity (or several nonlinearities 
conveniently located so that they can be incorporated in one de- 
scribing function) . These techniques can define stability and 
estimate relative stability for fairly complex systems as long 
as the conditions of nonlinearity are not too complex. Such 
cases areeasily treated using algebraic methods, the effect of 
the nonlinearity being represented as a movement of the operating 
point on the parameter plane, which in turn represents a varia- 
tion of the characteristic roots as a function of signal ampli- 
tude. The algebraic methods are capable of extending such analy- 
sis to systems containing two distinct nonlinear components, and 
can be used to predict the transient response of the system 
rather accurately. Techniques for such problems are developed 
m Chapter II. 

Chapter III is concerned with a much more difficult non- 
linear problem, that of asymmetrical nonlinear oscillations. 

These are oscillations consisting of a limit cycle superimposed 
oil another signal. The problems studied on the parameter plane 


involve steady— state operating conditions (rather than transient 
conditions) , and permit analysis of the existence of oscil- 
lations as well as their dependence on parameter values and in- 
put signal values. Extension to linearization with either 
signals is included, as well as some design considerations. 

It is felt that the results presented here indicate the 
capabilities of the algebraic methods in dealing with complex 
linear and nonlinear problems. It is also felt that the re- 
sults presented here will be directly applicable to a number of 
practical problems, and will point out avenues of approach to 
still additional problems. 
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SOLUTION OF EQUATIONS WITH COEFFICIENTS 
THAT ARE QUADRATIC IN « and B 

1.1 INTRODUCTION" 

It has been shown that the characteristic equation can be 
solved for a = <*();, «0 n ) and B = ^(C,^) when the coefficients of 
the characteristic equation are of the forms: 
a) aj. = b^a + c^H + d^, . 

h) ^ = V + °k* * V 6 + \ (!) 

c) °k = \ 2 a2 + b kl" + V S + c kl^ + C k2^ 2 + \ 

V + b k(n-l ) 0,11-1 + * • • \(n-l)“ n_iB + * ’ ' 

^(n-l)*"" 1 + c kn^ n + *k 

In addition practical solutions have' been obtained for the first 
two of these coefficient forms, i.e., computer programs have been 
written for them and successfully applied. The development to 
be presented here is a particular solution for case 1-lc, parti- 
cularly in the sense that a computer program has been obtained 
which solves the equations of a third order system for which the 
coefficients are quadratic m » and B, but which do not contain 
all of the or and B combinations indicated. At the same time the 
solution is a general solution in the sense that the -program can 
be modified to solve the equations of an n tb order system, and 
can also be modified to accept all of the o’ and B forms indi- 
cated m 

“k = b k2 a2 + b ki* + V s + C k!« + c k2* 2 + a* 

The modifications to be made in the program are discussed, h — 
the necessary programming has not been done. 


1.2 THE PROBLEM: Cascade Compensation with ‘two identical 

filter sections. 


In the design of feedback control systems it is common to 
use compensators which are filters placed in cascade with the 
main transmission path. Frequently two sections of filter are 
needed, and if identical sections are used with an isolation 
amplifier so that their transfer functions can be multiplied, 
then manipulation of the transfer function equation provides a 
characteristic equation m which" the coefficients are quadratic 
in z and p, the zero and pole of the compensators. For example 



(1-2) 

(1-3) 


14Q a =0 = 1+ . K(s +2zs+z j 

(s J +Xs'+Ys) (s*+2pS+p^) 


from which the characteristic equation is 

s 5 + (X+2p)s 4 + (p 2 +2Xp+Y)s 3 + (Xp 2 +2Yp + I<)s 2 + 

+ (Yp 2 + 2Kz)s + Kz 2 = 0 (1-5) 

A & 

Letting p = « and z * J it is noted that all of the forms speci- 
fied in the quadratic case definition of a^ do appear in at least 
some of the coefficients except that there is no o/B product term. 

The formulation gust given does not conform to normal con- 
trol system practice, however, in that an important restriction 
on the design of the compensator is the usual requirement that 
steady state accuracy must be maintained by keeping the error 
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coefficient unchanged. To do this the physical adjustment is to 
alter the gam of the amplifier, hut in the mathematical analysis 
it is more convenient to include this restriction in the transfer 
function of the compensator by defining (for this case) 


. <£\ 2 

K z J 's+p' 


‘This alters the algebraic form of the characteristic equation 


which becomes! 

0 = 1 + 


K(f) 2 (s 2 +2zs+z 2 ) 
(s 3 +Xs 2 +Ys) (s 2 +2ps+p 2 ) 


= (s 3 +Xs 2 +Ys) (s 2 +2ps+p 2 )+K S j{s 2 +2zs+z 2 ) 


= s 5 +(X+2p)s 4 + (p 2 f2Xp+Y) s 3 + [xp 2 +2Yp+K(E-) 2 ]s 2 

+[yp 2 +2Kp (|) ]s + Kp 2 (1-7) 

A „ ^ 

Choosing p = 8 and ~ this becomes 

0 = s 5 + (X+28) s 4 + (8 2 +2X8+Y)s 3 + (X<? 2 +2Yfl+Ka' 2 ) s 2 

+ (Y8 2 +2K8<*)s + K8 2 (1-8) 

In equation 1-8 the coefficients are quadratic in a and 8, 
but there is no term of the form b^®, and the program as written 
does not make provision for such a term, though modification of 
the program to include it is not difficult. The problem to be 
studied, then is that of a third order system compensated with 
two cascaded identical sections of filter, and with the added 
requirement that the error coefficient be maintained constant 
at a predetermined value. 
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1.3 DERIVATION Of THE GENERAL THIRD ORDER SYSTEM 
RELATIONSHIPS 

The general third order system is defined by the transfer 


function 


G(s) = 


(s+A) (s+B) (s+C) 


which is a Type Zero system, but which can be changed to Type 
1, 2, or 3 by setting one or more of the poles to zero. The 
compensator transfer function, including the gain multiplier 


which maintains the error coefficient is 

2 2 2 
. _ ,j 2 , 2 ,siz , 2 _ p fs +2zs+z ) 

G C - W - z 2 (s 2 +2ps+p2) 

From 1-9 and 1-10 the characteristic equation is 


( 1 - 10 ) 


[s 3 + (A+B+C)s 2 + ( AB+BC+AC ) s + ABc] (s 2 +2ps+p 2 ) + 


+ K ^(s^Szs+z 2 ) 
zr 


( 1 - 11 ) 


■This expands to 


S 5 + (A+B+C+2p)s 4 + [aB+BC+AC+2p(A+B+C)+P 2 ]s 3 

r 2 

+ [aBC+2p (AB+BC+CA) + p 2 (A+B+C) + K^js 2 + 

+ j^2pABC+p 2 (AB+BC+AC) + 2Kp (^) Js+p 2 (ABC+2K) - 0 


Let p - /3 


H A a 

z 


sum of root products talcen 2 at a time 


A+B+C ~ = sum of roots (poles) 

AB+BC+AC = I>i = sum of root products taken 2 at a time 

£ Hii - sum of root products taken n at a time 

ABC . . = 17 r^ - products of the roots 
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Then equation 1-12 becomes: 

s 5 + (£r.+ 2P)s 4 + (£ n r x + Z^P+P^s 3 + 

(H ^+2^ n r i P + ^}: i S 2 + K* 2 )s 2 + (1“13) 

(2D r P + Y n r,P 2 + 2K<*P)s + (n r.+2K)P 2 = 0 
1 u 2 1 

Collecting ‘like terms in a and P : 

a 2 (Ks 2 ) + «P(2Ks) + p 2 (s 3 + £r.s 2 + £ n j^s + II r ± +K) + 

+ P {2s 4 t 2^r i s 3 + 2^ n r^s 2 + 2^3) 

+ (s S + ^s 4 + £ Or^s 3 + Or.s 2 ) = 0 (1-14) 


Using the basic parameter plane relationships: • 
n 


1 <■-»* ^ = 0 



(1-15) 

"k=0 




n 




X (-D k ^ = O’ 



(1-16) 

k’=0 




and defining: 




B 21 = KW 2 ^!?) 



(1-17) 

B 22 = KU>2u 2 W 



d-18) 

D 2 = -2KWU 0 (C) 



(1-19) 

D 2 = -2K®U 1 (e) 



(1-20) 

E n = 2 u)4 U 3 (?) -2^0^02(5) 

+ 2 l 

DrA(E) - 
2 1 1 



- 2n 

ri «u 0 (0 

(1-21) 

E i2 - 2 " 4 ve> - 2 I r i“ 3u 3 ( ^ 

4 2 l 

H r.o: 2 U 2 (C) 


-20 r^U-^C) 



(1-22) 
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F 21 

= -« 3 u 2 (i;) + }_ - t 

l 2 1 r i UJU o (r ) 



+ ( n r jL +K)U_ 1 (0 


(1-23) 

F 22 

= -» 3 u 3 (a + X r ± “> 2 u a (r) -, 

'i n r^^tn 



+ (Or^+KjU^C) 


(1-24) 

G 1 

= -“ 5 u 4 (C) 4 lr.u) 4 u 3 (0 -I 

n r,w 3 U,(e) 



+ I^wV^S) 


(1-25) 

G 2 

= -n> S U 5 (C) + ^W i iu 4 U 4 (ir) -£ 

0 r.m 3 n (0 
2 i 3 



+0 r i <0 2 O 2 (C) 


(1-26) 

P 1 

= * D l 


(1-27) 

P 2 

= Pd 2 


(1-28) 

G 1 

= PE 11 4 p2p 21 4 G 1 


(1-29) 

°2 

= PE 12 4 ^22 4 G 2 


(1-30) 

This results in 




« 2 b 

21 4 °' P 1 4 Qi = 0 


(1-31) 

C 2 B 

22 4 ^2 4 Q 2 = 0 


(1-32) 


which are two non-linear algebraic equations completely general- 
ised in terms of the uncompensated system poles and root locus 
gain, «> and the first kind of Chebysbev Functions. These 
must be solved simultaneously for the correct values of a and 
P. To do this, the method with the best chance of success 
appears to be Sylvester's Method in which we form a set of four 
equations by taking the original Equations (1-31) and (1-32) and 


forming two more by a multiplication with o' giving: 



+ o-P^ 
+ ffp 2 


+ Q-L ■= 0 

4 a 2 “ 0 
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(1-33) 

(1-34) 



a- 3 B 21 + « 2 P 1 + 


<x 3 B 22 + « 2 P 2 + <VQ 2 = 0 


Now placing these equations in matrix form: 


0 

B 21 

F 1 



o' 3 

0 

B 22 

P 2 

°2 


a 2 

B 21 

P 1 

Q 1 

0 


a 

_ B 22 

P 2 

°2 

0__ 


1 

not 

zero 

then: 



0 

B 21 

P 1 

Q 1 


0 

b 22 

P 2 

q 2 


B 21 

P 1 

Q 1 

0 

= 

0 

B 22 

P 2 

«2 

0 



= 0 


(1-35) 

(1-36) 


(1-37) 


(1-38) 


Expanding this determinant 


" B 21 Q 2 4 B 21 B 22«1 Q 2. 4 . P 1 P 2 B 21 Q 2 


P 1 Q 2 B 22 


4 Q 1 Q 2 B 21 B 22 


-«1 B 22 - Ql P 2 B 21 + Q 1 P 1 P 2 B 22 = 0 < 1 " 39) 


Substituting equations (1-27) through (1-30) in equation (1-39) 
provides a fourth order equation in B : 

B (~ F 22 B 21 4 2F 21 F 22 B 21 B 22 4 D 1 D 2 F 22 B 21 ~ -F 22 D 1 B 22 


P 2 B 2 
F 21 B 22 


F 21 D 2 B 21 + D 1 D 2 F 21 B 22^ + 


B 3 (-2E 12 F 22 B 2 1 ^ 2E.,,F„,B, 1 ,B. 


ll i 22 I5 2l !: ’22 + 2P 21 E 12 B 21 B 22 4 


D 1 D 2 E 12 B 21 


D 1 E 12 B 22 


2E 11 F 21 B 22 ' " 


D 2 E 11 B 21 4 D 1 D 2 E 11 B 22 ) 4 

B (-E 12 B 21 “ 2F 22 G 2 B 21 4 2E 11 E 12 B 21 B 22 4 2F 21 G 2 B 21 B 22 


2G 1 F 22 B 21 B 22 4 D 1 D 2 G 2 B 21 ” D 1 G 2 B 22 " E 11 B 22 
2F 21 G 1 B 22 “ D 2 G 1 B 21 4 B 1 D 2 G 1 B 22 ) 4 


P <" 2E 12 G 2 B 21 4 2E 11 G 2 B 21 B 22 4 2G 1 E 12 B 21 B 22 " 2E 11 G 1 B 22> 4 


*~ G 2 B 21 4 2G 1 G 2 B 21 B 22 G 1 B 22* _ 0 


(1-40) 


from which the coefficients may be determined by a substitution 
of (1-17) through (1-26) and the values of the first hind of 
Chebyshev functions in terms of £ and w. since the solution of 
a fourth order equation is at best difficult, it is at this 
point a digital computer becomes a necessity. 

The major problem is not the actual solution of the quartic 
itself, but rather the proper choice of one of the four solutions. 
There are two marked characteristics, however, which help in the 
selection. These are: 

a) Complex answers to the quartic have no physical signifi- 
cance and may therefore be discarded as erroneous. 

b) The definition of a requires that a and 3 be of the 

same sign so that p and z will be of identical sign. 

, 2 
Using this information and that available from the Ross-Warren 

method as to compensator pole and zero location, it is found that 

the solution to the B quartic is the largest, positive, real 

value. 

Now entering equation (1-27) with this value, and evaluating 
the other coefficients 

» = [-Q 1 /B 2 i-l iS (1-41) 

for in the third order case P^ is always identically zero. 
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Thus, with the programming of -the appropriate equations, 
the digital computer could give all of the values and plot the 
constant zeta and constant omega loci on the parameter Plane for 
any desired values. 

1.4 SOME APPLICATIONS OF THE PROGRAM 

Several third order systems were investigated by the appli- 
cation of the generalized equations and the Parameter Plane 
curves. Figures 1-1 through 1-8 were plotted. Of these, the 
K/s 3 family appears the most interesting. Further investigation 
of three of the curves in this family, Figures 1-1, 1-2 and 1-3 
shows that there is a relationship between K, the root locus 
gam, a and P. 

These relations are: 

’ 3 

a) Choose a point on the 1/s e-P plane. 

b) Zeta reads directly. 

c) Determine the actual omega at that point by multiplying 
the value read by the cube root of the uncompensated 
system gain. 

d) Read the "a.tue of of directly from the point chosen. 

e) Read the value of P from the point chosen. 

f) Obtain the true value of P by multiplying this value 
by the cube root of the uncompensated system gain. 

By this method, the values of O' and P may be determined for 

V 

all — j systems from one universal curve, 
s 

1.5 BANDWIDTH CURVES ON THE Qf-p Plane 

In many instances, there is also a bandwidth criterion 

9 


imposed on the engineer as well as an optimal operating point for 
the plant under consideration. With this in mind, equations for 
the plotting of constant bandwidth curves on the tf-P plane are 
developed. For the purpose of this development a constant band- 
width curve will be defined as: 

A constant bandwidth curve for G(ga^) = M is a curve drawn 
upon the parameter plane which specifies the relation between 
the parameters necessary if the transfer function G(s), which 
is a function of the parameters, is to have magnitude M at 

5 

the real frequency 

Once these curves are obtained they may be superimposed on the 


parameter plane thus indicating what values of the parameters are 
necessary in order to meet the specifications. 


Taking the rational transfer" function and defining it: 

, w.i v" -w "" 1 * b! * 

’ QW - ■ , % , ... v ^ 


where the P m 's and q^s ate of the form: 

p u » g u “ 2 + v + v* + v + y 2 + 3 u 


U — 0/1^2,.. »fITl 

q v = a v » 2 + V + c v «P + d v P + e v P 2 + f y 



(1-44) 

(1-45) 
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Employing Equation 1-45 m the parameterized form the generalized 
compensated third order transfer function is: 

■w - Hfj- a-46) 


2 2 

= ot KS + 

2eBKs + B 2 K 



(1-47) 

= p 2 | 

IV + 

Y n 

2 , ' 
r . s + , 

) n r s + 


r.] + 

L 

4 - 1 

1 > 

- 2 1 

L 3 

1 -J 

B 

[ 2s 4 + 2 } 

n r s 3 
- 1 1 

+ 2 Y n 1 
A 2 

■V 

1 

+ 2 I ” r i s l + 


[s 5 + 

y n 

4 ' 

r.s + 

) n r.s 3 

+ 1 

n r s 2 l (1-48) 


L 1 

1 i 

"2 


3 1 J 


Making the definitions : 

. n 

A r = I etc - £ox W D r' E r' F r 

v=0 
even 
n 

\ = 'i <-D %(v “ 1)i ^ a v ; etc * for B i' c i' D i' E i' F i 

V=0 

odd 

m 

G r = I (-l)^ u u^g u ; etc. for WW 1 !: 
u=0 
even 

G x = I (-D 2{U_1)l0 b9u ; et0 * for H i' X i’ J i' K i' L i 


(1-49) 


and substituting in Equation 1-46, 


(<* G r 4 K r ) +3 (ggl^) ■ 

G(jl “ b) “ (p 2 D r + PE r + F r ) + + BE i + Fj 
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It follows that 

M 2 P^P 4 + m 2 q 1 ? 3 + (m 2 r x - r 2 )S 2 + m 2 v 1 P + 

(M 2 W 1 - W 2 ) =0 (1-65) 

Since the Parameter Plane for compensation purposes has already 
been determined it is now a matter of taking the computed a 
values and substituting them along with a constant value of 
omega and M into Equation (1-65) and then solving the P quartic. 
This has as its solution the largest, real and positive value of 
the four roots as before. 

1.5- EXTENSIONS TO HIGHER ORDER SYSTEMS _ 

Although the work presented to this point has been limited 
to third order systems and the program written for this specific 
case, investigation shows that generalized equations may be 
written which will allow the extension of the program to higher 

’ til 

ordered systems. It can be shown for a given n order system 
with no zeros to be compensated with two identical sections of 
cascade compensation, that the characteristic equation of the 
system may be generally written as: 

s n+2 + 2ps n+1 + p 2 s n + (z 2 s 2 + 2pzs + p 2 )K + 

j=n D=n 

k=l k=0 

2 P I <1 n sk + \ 1 . V r i> s * + 

k=n 3 ' ^r 1 3 

j“l 3 

k=2 

k=n 

)“ <£ n r i) s k = ° d-66) 

k=n+l 3 
1=1 
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where for n=4 the equation would be written: 


+ 2ps 5 

+ p 2 S 4 + 

(z 2 s 2 + 2pzs + 

P 2 )K + 


2 p(Z 

n r.s 4 + 
1 1 

y n r s 3 + y n 

L 2 1 L 3 

r^s 2 + 

1 % r i s) + 

* 2 <r 

n r.s 3 + 
1 1 

y n r . s 2 + y n 

L 2 1 u 3 

r i s + ] 

> Ur.) + 
- 4 1 

6 n 

r.s 5 + T 

n r s 4 + y n r 

s 3 + y 

H r.s 2 ) = 

L x 

i Li 

2 1 ^3 1 

4 X 


(1-67) 

It may be further shown that the parameters defined by Equations 
(1-17) through (1-26) may be written: 


B 21 

= KwV^S) 

(1-68) 

B 22 

= KU> 2 U 2 (C) 

(1-69)' 

D 1 

= -2KWU 0 (?) 

(1-70) 

°2 

= -2KWU,($) 

3~n 

k=l 

(1-71) 

2(-l) n+1 w n+1 U n (4) + 2 Z [ (- 1 ) k " ku k_l(^ (I 11 r i' J 

(1-72) 


k=n 

j=l 


j=n 

k=i 

e 12 = 2(-l). n+1 w n+1 U n+;L (4) + 2 £ [(-l) k !^U k U) ( £ n r^] '(1-73) 

k=n 3 

j=l 

j«n 

k=0 

p 2i = I-dVu^io ■» [<-i) ls u? t u k _ 1 (S)(Z n r.) j + ku^U) 
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k=0 

P 22 = (-l) n «Pu n (S) + I [(-D kJC U k (€) (I r *i)] + (1 “ 75) 

k=n-l 3 

j=l 

k-2 

G x = (-l) n+ 2 ^ + 2 u n+1 (« + l [c-D^Vl^xI ^! 1 ] a ~ 76) 

k=n+l 3 

j=l 

j=n 

k=2 

Gs ■= (-l) n+2 ^2 Un+2(4) + l [(.D^teH^r^ (1-77) 

k=n+l 1 

3 - 1 

These then are the recursive equations required for the complete 
generalization to a n th order system. By employing the above 
equations and replacing in PROGRAM PROJECT cards 100 through 150 
and 300 through 540 with the appropriate programming, the program 
may be used for any given n* 31 order system. 

In like manner by generally defining: 


p (s) = a 2 KS 2 + a?Ks + P 2 K 


0 (s) = s n+2 + 2j3s lvfl + P 2 s n + 2 £ (I 17 X x )* + 

"L- =rn 


j-n D= n 

(1 ' 7 

k=n-l 3 k=n+l 3 

D”1 j =1 

and using Equations (1-49) through (1-52) we may replace in the 


program cards 2860 and 2880 through 2920, thus adapting this 
part of the program to a general n^ order application. 
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1.6 COMMENTS 

Throughout this development of the Parameter Plane quadratic 

extension, the c- 's in the generalized coefficient form; 
n K 

£ (h^u 2 + c^u + c^orp + e^P + ^P 2 + g k ) =0 (1-80) 

k =0 

have been identically zero. This at first appearance might seem 
to detract from the generalization. The inclusion of this para- 
meter does not however introduce any great difficulty in the 
solution. The change in the development would be to the value 
of and P 2 which would become: 

P 1 = C 1 + PD 1 (1-B1) 


P 2 = C 2 + 6D 2 (1-82 

and the final solution for o/ which would change to: 

r 2 -i % 

^1 ^1 " 2 ^ 21^1 

“ = - r - r - + — , -- (1-83 

"21 4B zl 

For this case, new selection rules for acceptable values of a 
would be used, and would be much like those presented for P. 

Though the extension of the Parameter Plane to include the 
n - 3 quadratic case makes this tool even more useful, further 
work is still to be done in this field. Not only must the 
equations for the solutions of the Parameter Plane curves for 


such cases as: 


a k = \" 2 e 2 + -l d^p 2 + e ^ 2 + f ^ 2 + 


g k »e + v + V + % 
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APPENDIX I 


PROGRAM PROJECT is designed to solve the or quadratic and S 
quartic. The program is divided into two main sections, the first 
for the computation of the i— P points and the second for the band- 
width points. 

The first section computes an 80 by 80 matrix of the and 
S points corresponding to set values of £ and w. The computational 
part is followed by two distinct graphing sections, one for lag 
and the other for lead compensation. 

The lag graphing section is set up so that during the plotting 
of the curves each value of or is tested to determine if its value 
is 10~ 7 s a s 1.0001. If no points are found within this range 
then a print out is made: 

NO LAG COMPENSATION POSSIBLE 

For the lead section graphs, or is again tested by the criter- 
ion 1.0001' s a u (X-graph scale) (X graph width). Again if there 
are no values of or within this region the statement: 

NO LEAD COMPENSATION POSSIBLE 

is printed. In this case however a study of the printed values 

of or must be made to insure that the points are indeed non-existant 

or rather just lie outside the range of the graph. 

* 

The second main section of the program computes the value of 
0 for a given value of or is determined by tbe X graph scale. Here 
the plotting routine is set up so as to not plot zero .points and 
to stop the curve when either the or or B value exceeds the range 
of the graph. 
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PROGRAM PROJECT 000000 

THIS PROGRAM COMPUTES THE VALUCS OF BETAtPOl.E LOCATION) AMD ALFAfPOLE 
—ZERO RATIO) BASED ON PARAMETER PLANE TECHNIQUES, THE COMPUTED , 

ALFA AND BETA VALUES ARE THOSE REQUIRED TO PLACE THE ROOTS OF ANY 



C IF IABCHP=1 CARDS FOUR THROUGH THIRTEEN ARE OMITTED 

C 

C«*«CAR0 FOUR - ILGPLT - COLUMN ONE IN FIXED, POINT 
C 0 - THE LAG ZONE CURVES WILL BE PLOTTED 



C i - THE LAG ZONE CURVES WILL MOT BE PLOTTED 

C ****** 

C IF IL6PLT--1 1 1 IE NEXT FOUR CARDS ARE OMITTED t 
C- ****** ' - 

C***CARD FIVE -IT(1)-ITI6) - COLUMNS 1-40 IN -ALFANUMERIC CHARACTERS' 

C this IS THE FIRST LINE OF THE LAG GRAPH TITLE 

Ci*<CARD SIX - IT 1 7 1 -I T 1 12 ) - COLUMNS 1-40 IN ALFAMUMERIC CHARACTERS 
C THIS -IS THE SECOND LINE OF THE LAG GRAPH TITLE, 

C***CARD SEVEN - L8L1 M J-LBLI20 ) - FOUR COLUMNS PER LABEL (TEN LABELS- IN 
C CONSECUTIVE COLUMNS) IN ALFANUMERIC CHARACTERS, 

C ' THESE ARE THE LABELS TO BE PUT ON THE CONSTANT OMCGA CURVES. TO 
C DETERMINE WHICH VALUES WILL BE PLOTTED, DIVIDE WFIN BY. 10 . THIS 

C ' .VALUE AND INTEGER MULTIPLES OF IT TO 10 WILL BE PLOTTED. 

C**»CARD EIGHT - XLGZ » YLGZ ' - TEN COLUMNS PER NUMBER IN EXPONENTIAL OR 
C FLOATING POINT. 

C THESf ARE THE X AND Y SCALES FOR THE' LAG- GRAPH. ONLY ONE SIGNl- 

C FICANT NUMBER IS TO 8E USED. 

C***CARD NINE - ILDPLT - COLUMN ONE IN FIXED POINT 
C 0 - THE- LEAD CURVES WIL'L BE PLOTTED 

C 1 - THE LEAD CURVES WILL NOT BE PLOTTED 

C ****** 

C IF ILDPLT=1 THE NEXT FOUR CARDS ARE OMITTED 

C ****** 

C* **CARD T.EN - THE SAME AS CARD FIVE EXCEPT FOR THE LEAD GRAPH 
C***CARD ELEVEN - THE SAME AS CARD SIX EXCEPT FOR THE LEAD GRAPH 
C* **CARD TWELVE - THE SAME AS CARD EIGHT EXCEPT FOR THE LEAD GRAPH 
C** "CARD THIRTEEN - A**<DUPLICATE*** OF CARD SEVEN 
C***CARD FOURTEEN - IBWCMP - COLUMN ONE FIXED- POINT 
C O' - BANDVIIDTH COMPUTATIONS AND GRAPHING WILL NOT BE DONE.. ■ 

C 1 - BANDWIDTH COMPUTATIONS WILL BE DONE ' 

C ****** 

C IF I BWCMP -0 THE REMAINING CARDS ARE OMITTED 

C ****** . • 

C***CARD FIFTEEN — 'BWXtBUY - THE SAMEAS CARD EIGHT - EXCEPT FOR, THE 
C ' BANDWIDTH CURVES. ' - 

C BWY IS ALSO USED. TO DETERMINE WHICH VALUES OF ALFA WILL BE USED IN 


c THE BANDWIDTH COMPUTATIONS. . 

C***CARD SIXTEEN - WEND - TEN COLUMNS IN FLOATING POINT 
C THIS IS THE MAXIMUM VALUE OF OMEGA FOR WHICH THE BANDWIDTH 

C COMMUTATIONS WILL BE DONE 

C* X* CARD SEVENTEEN - IBWPLT - COLUMN ONE IN FIXED POINT 
O'- THE BANDWIDTH CURVES WILL BE PLOTTED 
1 - THE BANDWIDTH CURVES WILL NOT BE PLOTTED 
****** 

IF I0WPLT=i THE REMAINING CARDS ARE OMITTED 

■K , 

C*-**CARD EIGHTECN - THE SAME AS CARD FIVE EXCEPT FOR THE BANDWIDTH CURVES 
C** *CARD NINETEEN - THE SAME AS CARD SIX EXCEPT FOR THE BANDWIDTH CURVES 
C*** CARD TWENTY - BANDWIDTH CURVE LABELS 

TO DETERMINE WHICH CURVES WILL BE PLOTTED, DIVIDE WEND BY 10. . 

THE PROGRAM PLOTS THIS CURVE AND INTEGER MULTIPLES OF IT UP TO 10. 


C 

c 

c 

c 

c- 


IT IS RECOMMENDED THAT FOR THE INITIAL RUN THE FOLLOWING DATA CARDS 

. BE CARDS 1 , 2 , 3 ( I ABCMP = 0 ) » 4t ILGPLT=1) ,9 ( ILDPLT = 1) , 14( IBWCMP =0 ) 

THESE DATA CARDS WILL ALLOW ONLY THE ALFA-BETA COMPUTATIONS TO BE 

COMPLETED. A PRINT OUT OF THE VALUES WILL BE OUTPUT WHICH .WILL ALLOW 
YOU TO CHOOSE THE PROPER CURVES AND SCALES. CAREFUL SELECTION 
OF CURVC SCALE5 IS IMPORTANT, FOR THE PROGRAM WILL NOT' ALLOW POINTS 
OUTSIDE THE AXIX LIMITS TO BE PLOTTED. 


DIMENSION AFI1K 00,00) ,BFIN ( 80, 00 1 ,XAZ ! CO ) >YBZ(BQ) sXAW(OO) , 

1 YBW180) ilTUZI ,LBL(20) .BCOri [ 5 ) ,ROOTR 1 4 ) ,ROOTl I/ t ) ,ACOFI(3) , 

2 UI10) .AROOTI (4),ACOFR<3) , BCOFR 1 5 ) ,ULABI 00 ) »ZLAB 100) ,AR00TRS4) 
COMMON BCOFR, BCOFI a ROOTR«ROOT I , BFINAL* ) FLAG, AF IN,BFI H 

9999 PRINT 140 
140 FORMAT tlHl)' 

DO 60 JK=1 » 6400 
AFINtJK) =0.0 
60 BFIN(JK) = 0.0 


000010 

000020 

000030 

000040 

000050 

000060 

000070 

000000 

000090 



READ 1 , A > B i C > G 

1 t OllHAT (4F10*0) 

TROD = A"B*C 
"UM = AiB+C 

•PRO “ A<B + A*C + B«C 
i-UDGM = PROD +■ G 
ZETA = 0.0 
READ 2.WF1N 

2 FORMAT (FIO.O) 

READ 9. IABCMP 

9 FORMAT (tl) 

IF(IABCMP-1>23,24,24 
23 STP = WFIN/80.. 

DO 12 L = 1,30 
t_J = SOML-l) 

V! = STP 
30 Utl)=-1. 

U( 2 ) =0, 

U(3>=1. 

DO 10 N=2,6 

10 U(N+2)=2.*ZETAX-U(N+1>-U(fi) 
DO 11 J=l,80 
LJ = LJ + 1 
Y/2=\/*W 
U3=W2*W 
W4=!/2»W2 
W5=W2*W3 
CORN = G*W2 
CON = -2.*G*W 
CONI = 2. *1/4 
C0N2 = ~2.*SUH*U3 
COM3 = 2.*ShPRD*W2 
COM4 = -2.*PR0D*ll 
C0N5 = SUM* W2 
C0N6 = -S/1PRD*W 
COM 7 = S'Ji 1* 1/4 


000100 
00 " 1 10 
000120 
000130 
000140 . 

000150 
000160 
000170 
000180 
000182 
000185 
•000188 
000190 
000200 
000210 
000220 
000230 
000240 
000250 
000260 
000270 
000280 
000290 
000300 
000310 
000320 
000330 
000340 
000350 
000360 
000370 
000380 
000390 
000400 
000410 
000420 


CONS = ".PRDSl/3 
CO/19 = PROD*U2 
B21 = COK!!'Ut3> 

B22 = CONN-U( 4) 

D1 = CON"U( 2 ) 

D2 = C0K*l!<3) ' 

£11 tz C0N1 J -U( 5 ) + C0N2*U14) + CON3 J U(3) + CON4*U(2> 
E12 = C0!11*U16) + CON2*U/5) + CON3*U(4) + CON4*U(3> 
F21 = ~U3*U(4) + CONS'U 1 3 > + CON6*U(2) + PRDGN*U(1) 
F22 = -’J3*U(5) + CON 5* U 1 4 ) I C0N6*U(3) \ PRDGH’U<2> 


61 = -U 5 s U ( 6 ) + CON7*U(5) + COK8~U(4> + 

62 = -W 5*0(7) + CON7*U(6) + CON8*U(5> + 


CON9‘”U 1 3 ) 

C0N9*U(4) 

COF1 ="B21*F22*(2.*F21>'822-F22*B2l)--F21* ( F21*022"B22-! 02*02*021 ) 
C0F2 = Ell* (2«*B22*(F22-' t B21-F21 I B22)-D2<D2 :; B21H-2.*E12*B21* (F21*82 
12-F22*B21) T ' 

COF3 = B21* {-821* t E12-E12-S-2.* F22 *G21 -D2*D2*61+2 »*B22* t Ell'£12+F21* 
1 G2-/G12F22) )-B22*B22*<Ell*EU+2»*F21<Gl) 

COF4=2.*G2*B21* S Ell* B22-E 12*321 )•< 2 .*B22*G1* { El 2* 8 2 1-El 1*022 > 

COF5 = -(G2*B21-G1*822)*1G2*B21-G1*B22 > 

DO 50 I =1,5 
50 BCOFIU) = 0.0 
BCOFRdl = 1.0 
BCOFR1 2 > = COF2/COF1 
BCOFR 1 3 1 = C0F3/C0F1 
BC0FR{ 4 1 = COF4/COF1 
ocorR(5) = cors/coFi 
CAUL ABCTART 
I FLAG = 0 
CALL SORT 

IF I IFLAG-1 1300,11,11 
300 BFtfULJ) = BF INAL 

Q1 = GF IIH LJ) * ( Ell+OFIN I LJ 1* F21 i+Gl 

ACOFRI 11=1.0 

ACOFR1 2 >=0.0 

ACOFRt 3 ) = Q1/B21 

ALFASQ = ABSF ( ACOFR (3)1 


000*30 
000440 
000450 
000460 
000470 
000430 
000490 
000500 
000510 
000520 
000530 
000540 
000550 
000560 
000570 
000580 
000590 
000600 
000610 
000620 
000630 
000640 
000650 
000660 
000670 
000680 
000690 
000700 
000710 
000720 
000730 
0007 ’.0 
0<’.J?30 
C. >760 
0.10770 
000780 



AFIN(LJ) - SORTF ( ALF ASO) 000790 

U U = W+STP 000000 

12 ZETA = ZETA + .0125 000010 

L3l.fl) = 4)SZ«*0 000020 

L3L(2) = 4liZ--.il 000330 

' LBLCD = 4liZ = *2 000840 

L8L (4) = 4HZ- .3 000050 

LBL ( 5 ) = 4HZ=»4 000860 

LBL ( 6 ) 4!IZ=.5 000870 

LBL (7) 4HZ=*6 000880 

LBL (8) = 4HZ-* 7 000890 

LBL<9] = 4HZ=«8 000900 

LBL (10) = 4HZ=.9 000910 

READ 7, ILGPLT 000920 

7 FORMAT (11) 000930 

IF{ ILGPLT-1)8>67>67 000940 

8 READ 3,- UTt t)»I = l>12' 000950 

3 FORMAT (6A0) 000960 

READ 6* ( LBL < I ) » 1=11 i20) 000970 

6 FORMAT (10A4) 000900 

READ 4. XLGZ * YLGZ 000990 

4 FORMAT <2E10„0) ,001000 

XLGLM = 9,*XLGZ 001010 

YLGLM = 15«*YLGZ 001020 

MODE = 1 001030 

IL = 0 001040 

DO 62 K=1»C0>8 001050 

LL = 1 001060 

KJ = (K-11500 001070 

DO 61 J=1 9 80 001080 

KJ = KJ+1 001090 

IF<AFt(UKJ)-*0000001 161,6110 >61 10- 001095 

6110 IF(AFIN(KJ)-1, 0001)6113)61, 61 001100 

6113 IFIAFIIKKJ) - XLGLM) 6114*61 j61 001110 

C CARDS 1120 - 1130 ARE MISSING 

6114 XAZ(LL) = AFINtKJ) 001140 


I F < BF I N ( KJ ) - YLGLM) 61 12*61*61 001150 

C CAROS 1160 - 11.0 ARE MISSING 

6112 YBZ(LL) = BrtN(KJ) 001100 

LL = LL + 1 001190 

61 CONTINUE 001200 

JJ = LL - 1 OOT210 

IL = IL + 1 001220 

I F ( JJ-1 )62* 62 .6116 001230 

6116 LAL = LBL(IL) 001240 

CALL DRAL'f JJ *XAZ s YBZ *MODE * 0* LAL > IT >XLGZ* YLGZ *0* 0* 0 *0*9*15*0* LAST ) 001250 

6111 MODE = 2 001260 

62 CONTINUE 001270 

IFUtODOl 165*65*6120 001280 

6120 DO 66 K= 8 1 00 *8 001290 

LL = 1 001300 

DO 63 J-l * 80 001310 

JK = (j-lKOO + K 001320 

IF (AFIN(JK) -*0000001 163,6127,6127 001325 

6127 IFIAFliKJKl-l. 0001)6123. 63|63 001330 

6123 IF(AFINIJK) - XLGLM) 6124 * 63 * 63 001340 

C CARDS 1350 - 1360 ARE MISSING 

6124 XAU(LL) = AFIN(JK) 001370 

IF< BFI N( JK ) - YLGLM) 6122 *63* 63 001300 

C CARDS 1390 - 1400 ARE MISSING 

6122 YBW(LL) = BF IN ( JK) ' 001410 

LL = LL + 1 001420 

63 CONTINUE 001430 

JJ = LL - 1 001440 

IL = IL + 1 001430 

IT! JJ-1)6121»6121,6126 001460 

6121 I F ( K-30 )66 ,6125,6125 001470 

6125 MODE = 3 ' 001480 

. L/L = 4H 001490 

JJ = 2 001500 

XAUU) = XLGLM 001510 

XAWI2) = XLGLM 001520 



YGL'El) =0.0 
YBUE?) = Yl.GZ 
GO TO 2000 
6126 LAL = LBLEIL)' 

2000 CALL DR JJtXAW, YOU , MORE .0 . LAL* IT >Xl GZ» YLGZ ,0,0, 0,0, 9, 15,0, LAST ) 

• MODE =2 •• 

IKE £-■ i '66,64,64 

6 A MODE = . 

66 CONTINUE ' . . 

• GO TO 67 
65 PRINT 130 

130 FORMAT UX.33H NO LAG COMPENSATION IS POSSIBLE »///) 

67 READ 20 » ILDPLT ‘ ■ 

20 FORMAT (III 
IFtILDPLT-1160,1000.1000 

63 READ 5, ( I T ( I ) , 1=1 , 12 ) 

5 FORMAT ( 6A6 ) 

■ READ 21. XLDZ.YLDZ 

21 FORMAT (-2E10.0) 

READ 22, (LBL< II ,1=11,20) 

22 FORMAT (10A4) 

XLDLM = 9.<-XLDZ 
YLDLM = 15.S-YLOZ 
IL = 0 

MODE = 1 

DO 72 K = 1,80,8 
KJ = (K-l)*80 
KK = 1 

DO 71 J » 1,80 
KJ = K.J + 1 

IF t AFINOCJ) -1 .0001 17 1,7 ill ,7111 
7111 IF ( AFI/I E KJ1 - XLDLM)7117,71,7l 
,C CARDS 1830 - 1040 ARE MISSING 

7117 XAZEKK) = AFINEKJ) 

IFEBKINEKJ) - YLDLM) 7118,71,71 
C CARDS 1870 - 1880 ARE MISSING 


001530 

001540 

001550 

001360 

001570 

.001580 

001590 

001600 

001610 

001620 

001630 

001640 

001650 

001660 

001670 

001680 

001690 

001700 

001710 

001713 

001716 

001720 

001730 

001740 

001750 

001760 

001770 

001700 

001790 

001800 

001810 

001820 

001350 

001060 


7110 YiV.EKKI = BFKIEKJ) 001890 

KK = KK + 1 001900 

71' CONTINUE 001910 

KM = KX-1 001920 

IL ' IL H 001930 

1FEKM-1 172,72,7119 001940 

7119 LAL = LEM. ( IL) - 001950 

CALL DRANEMMsXaZ , YBZ ,MODC , 0 ,LAL , IT , XLDZ , YE.DZ ,0,0, 0,0, 9, 15,0, LAST ) 001960 

7110 MODE = 2 001970 

72 CONTINUE 001900 

IF (KODE--1 1.70, 70,78 001990. 

78 DO 76 K=8,80,8 YO2000 

KK = 1 002010 

DO 73 J = 1,80 002020 

JK = (J-l)*80 + K 002030 

IF ( AFt M ( JK)— 1.0001 (73,7121 ,7121 002040 

7121 IFEAFIN(JK) - XLDLM ) 7,127, 73 ,73 002050 

C CARDS 2060 - 2070 ARE MISSING 

7127 XAl/EKK) = AF1NEJK1 ' . 002080 

IFIBFINEJK) - YLDLM) 7128 ,73,73 002090 

C CARDS 2100 - 2110 ARE MISSING ' * " 

7128 YBWEKK.) = BFINtJRl 002120 ' 

KfC = KK + 1 ■ 002130 

73 CONTINUE 002140 

MM = KK-1 002130 

IL = IL'+ 1 002160 

II-IMM-l (7120,7120,7) 29 002170 

7120 1KEK-80176, 7122, 7122 ' 002100 

7122 MODE = 3 * 002190 

LAL = 4H * 002200 

Mil = 2 002210 

XAUEIJ = XLDLM 002220 

XAUE2) = XLDLM 002230 

YBW(l) = 0.0 002240 

YE5UE2) = YLD7. 002250 

GO TO 2001 002260 



73 

2 00 3 


75 

76 

70 

131 

1000 


Cl 

02 

100 

101 

102 

103 

in 

112 

333 

’4 

314 

210 

219 


LAI = Lni.t il 1 . , _ 

C,M 1. pRAU<n.i,XAl'<YBU»M0!'E»0,l AL,IT,Xl i>. 

Mul'T “ 2 

1 1" ( K— 72 1 7 6 , 7 5 ,75 
mode = 3 
torn inuc 

1-3 TO 1000 
( Hit 131 
i ORMAT (IX, '.H NO 
CONTINUE 
ZLAB(l) = 0.0 
DO 81 3=2,10 
ZLAB1 * 1 = ZLABII 
ULAB ( 8 ) -• Oc*ST!> 

DO C2 N-16»80s8 
ULAG(il) = ULABIK-8) + WLAB(B) 

Pi; HIT 100 
FORMAT (1111) 

PRINT 101 

FORMAT ( 2X, 20H THE ALFA VALUES ARE 
PRINT 102 * tZLA..<I! 1=1,10] 
rORMAT 11X,6H ZETA '0F11.6) 

PRINT 103 , CULABt J 1 , ( . MU J* 1 1 » I =1 >73* 8 ) J-8i 00 , 8) 
FORMAT (/»1X,F6,2,10U i ,5) 

FORMAT 1 t/////?X»20!l THE BETA VALUES ARE .//l 

PRINT 112, (ZLABt!) 1=1.101 
FORMAT (IX, 61! ZETA »10F11,6) 

PRINT 113. (ULABI J 1 , (BFIN ( J, I ) , I s ‘l »73» 8) J=8, 00 ,8 ) 
FOPvMAT ( / ,lA,r6«2. 10E1 1,5) 

PRINT 114 
FORMAT (1H1I 
READ 218, 1BUCMP 
FORMAT (II) 

IFt IBL'CMP-1 11002 ,219,219 
READ 217, BWXtBUY 


Yl D7 , 0 , 0 1 0 1 0 v 9 s 1 r > ,0 , LAST ) 


LEAD COMPENSATION IS PO5SI0LE ,///) 


- 1 ) + . 1 . 


,//) 


002270 

002233 

002.290 

002300 

062310 

002320 

002330 

002340 

002.390 

002300 

002370 

002300 

002390 

002400 

002410 

002420 

002430 

002440 

002450 

002400 

002470 

002480 

002490 

002500 

002510 

002520 

002530 

002540 

002550 

002560 

002570 

0025G0 

002590 

002600 

002610 

002620 


217 FORMAT 12210,0) 

READ 223, RENO 
223 FORMAT (F10.0) 

TBULM = 15.* GUY 
XLM = 9 ,*B1)X 
AM2=«5 

STEP = UEND/20. 

W = STEP 
ALFASP = XLM/20 « 

XAU ( 1 ) - ALFASP 
DO 200 K=2, 20 

200 XAU ( X 1 = XAU(K-l) + XAU(l) 

DO '’03 f!=l, 20 ‘ 

AL A = ALFASP 
DO 210 11=1,20 
XAZ(M) = 0,0 
210 YBZ(M) = 0,0 
DO 207 1=1,20 
YBU(N) = U 
K2=11*W 
113=11*1/2 
114=112*1(2 
115=112*113 
. ACR = -1/2*6 ’ 

- AXR = G. 

ADR=2.* 114-2 , #U2*SMPRD 
All = 2, *6*1/ 

AER=-1'2*SUM + PROD 

APR =114 -SUM - 112-fROD 

ADI=-2,>-U3*SUM * 2,M!*-PROD 

AC I = -1/3 + !/*si:p RD 

Arl = 115 - U-KSflPIlD 

PI = ADR*ADR:ADI*ADI 

01 = 2,*ADR*AER -f 2,*ACJ7ADI 

R1 = .2«*ADR*AFR + AER*AER + AEI*AE1 3 

VI = 2,*AER*AFR + ?,*AEI*AFI 


002630 
002633 
002636 
002640 
002650 
002660 
002670 
002680 
002690 
002700 
002710 
002720 
00273' 
0027 . 
002750 
002760 
002770 
002700 
002790 
002800 
002010 
002020 
002830 
002840 
0Q2850 
002060 
002070 
002000 
002090 
002900 
002910 
002920 
002930 
002940 

2«*-ADI*AFI . 002930 

002960 



1/1 = ArRi'ACR + AFI<Ar.I 
A2 = ALFAKALTA 
A4 =■ A? = A2 
R2 = AZ^AI I 'A 1 1 

1.1/2 = AAJ"AGR"-AGR +- 2.*A2£AKR«/iGR + AKR^AKR - 
DO- 51 11 = 1 ,5 
51 GCGFI (ft) = 0.0 
BCOFRtl) = 1.0 
3C0FR t 2 ) = Ql'/Pl 
BC0FR(3) = ( AM2FR1--R2 1 / { AI12-P1 ) 

BCOFRt'A) = VI /PI' 

6C0FR 1 5 ) = (AM2*Vfl-WW2)/(AH2*Pl) 

CALL ABETART 
I FLAG = 0 ■ 

CALI. SORT -.J 

IF ( 1FLAG-1 ) 900)206)206 . 

206 F.FIUtN.tl = 0.0 
GO TO 207 

900 BF1H1N.I1 a BFSNAL 

207 ALFA = ALFA + ALFASP 
‘ 203 U = M + STEP 

READ 216. .I3UPLT 
216 F.0R11AT till 

• IF 1 1B1/PLT- 11216 « 1001 . 101 
214 flODE = 1- 

READ 202 >(.IT(K)»K=1»12) 

202 TORilAT (6A8 1 

READ 201. ( LBL ( N) . N=2 > 20 >2 ) 

201 FORMAT (10A41 
DO 211 K-“2 .20.2 
KK = 1 

IF(N-20)2Q4. 205.20! 

205 MODE = 3 
204 CONTINUE 

DO 212 1=1.20 

irtSFUKN.n - .000001 1212.212.209 


002970 
002900 
.002990- 
003000 
003010 
003020 
003030 
003040 
003050 
003060 
00,3070.' 
003000 
.003090 
003100 
003110 
003120 1 
003130 
003140 
003150 
003160 
003170 
003100 
003190 
003200 
003210 
003720 
003230 
003240 
003250 
003260 
003270 
’003200 
CO .790 
003300 
003310 
003320 


209 IFtDrTNtH.il - YBULHJ 905 .906.906 

905 .yozocd a ormiN.i 1 

X’UOCK) « XAUt I ) 

CARD 3360 I c 
= K.K. + 1 


212 CONTINUE 
906 JJ = KK. - 1 

IFl JJ-1 1221 .221.220 

221 1 T C H— 20 1211.225.225 
225 IFU.0DE-D224.224.215 
215 MODE. = 3 

LAL =’ 4H 
XAZtll =0.0 
XAZI21 = 0.0 • 

YBZI11 = 0.0 
YBZtZ' = BUY 
JJ = 2 
GO TO 2002 

220 LAL = LBUH1 ' 

2002 CALL DRAW ( JJ .XAZ. YBZ .MODE .0 .UAL. I T s 8UX 

222 IF IN— 20 1208.211 .211 
-33 MOOS = 2 . 

211 CONTI HUE 
60 TO 1001 


.BUY. 0.0.0 >0.9. 15.0 .LAST 1 


226 FORMAT 2 Dill .1X»76HTH' : .RE ARE NO POSSIBLE BANDWIDTH CURVES FOR THE F 
’ 1IKAL VALUE OF OMEGA ST AT FD ' 

GO TO 1002 

FORMAT^ (1H1 .20X.38H VALUES OF BETA FOR CONSTANT BANDWIDTH 
PRINT 122 . (YOUtiD .K--2.20.2) 

FORMA 1 I/.9X.10F11.6) 

PRIM! 121. t\ ViltX) .tor ’.‘ft J.Kl . J=2.20.2) .K=l»20) 

FORMAT t/.l- *6.2.2X91 Ot: 11. 6) 

CONTINUE 
GO TO 9999 
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122 

121 

002 


003330 

003340' 

003350 

003370 

003300 

003390 

003400 

003410 

003415 

003420 

003430 

003440 

003450 ■ 

003460 

003470 

003480 

003490 

003500 

003510 

003520 

003530 

003540 

003541 

003542 

003543 

003544 

003545 

003550 

003560 

003570 

003580 

003590 

003600 

003610 

003615 



DI HEMS I of a' ( 5 1 Iy IHAG ( 5 1 ,U < 4) .V ( 4) , !l t 50 ) , B( 50,1 ,C t 50 ) , D < 50) ,Et 50) 
1 ,C<MV(50) 

Dlf’.EI.SIOIl AFIN< GO, SO) ,BFIM(S0s80) . 

CO! .1 iOIl A » Y It IAG 5 U * V » DUffMY 1 s DUMf ) Y2 , AF If! . 3F I H 

n - 4 

r= 10,0 

L-29 

IER--0 

• irttl! 54,54,52 
34 IER=i. 

52 KP3»tU3 
100 3(2)=0<,0 
B{ 1 1 =0o 0 
C(2)"-0.0 
c(i:=o,o 

D ( 2 : -o C 0 

E[2J“0„0 

K(2 J “Oc 0 

CO iOl J=3,IIP3 
101 H(J)=A(J“2) 
t=j «o 

SK=10oO"-"F 

150 1 F ( I! I RP3 1 1 200,151 ,200 
131 U<KP3)=0,0 
V(KP31 -OaO 

cc::v(.'5P3» =sk 
! p3- t:: j 3“i 

IF ( f Is' 3 ) 132-1 52 >150 
152 IER=1 

ZOO Ii r [KP3-3)205j51,201 
205 ICR=1 
201 P5=0o0 


003020 

00363 

00360 

003650 

003660 

003670 

003600 

003690 

003700 

00003710 

00003720 

00003730 

00003740 

00003750 

00003760 

00003770 

00203730 

00003790 

00003000 

00003810 

0000302.0 

00003830 

00003040 

00003850 

00003860 

00003070 

00003880 

00003390 

00003900 

00003910 

00003920 

00003930 

00003940 

00003950 

00003960 

00003970 


os-o.o 
PT--0o0 
01=0 „0 
5=0,0 
REY=1,0 
5K=10<,0<,F 
IF( RP3— 4) 206 s 202 »?03 

206 IER=1 

202 R=-ST('')/!H3) 

GO TO 500 

203 DO 207 J=3,t-!P3 

204 5 207 s204 

204 s-siLOGriAc-srintJ) n 

207 COHTiJiUE 
FP;!l=fl-H 
S-E-SPFtS/FPril » 

CO 208 J=3»f!P3 

208 H(J)=HU)/S 

210 IF( ABSFtHl 4 ) / !!( 3 ) ) —AG5F flit NP3-1 ) /H 1 MP3 )))250s252s252 

250 T=— T 
I1~-(NP3-4>/2 * 3 
DO 251 J=3»M 
C=HiJ) 

JJ=[!P3-J+3 
il(J)=H( JJ) 

251 iHJJ)=S 

252 I F ! 05 > 253, 25 'n 253 

253 P=PS 
Q--G5 

G 3 . 3 300 

254 Hl!2=ll( MP3-2 ) 

1F!!*!I2) 256,255,256 

255 0=1,0 
P=— 2,0 

GO TO 257 
.56 0“IHMP3)/HH2 


00003930 
00003990 
00004000 
00004010 
00004020 
00004030 
00004040 
00004050 
00004060 
00004070 
00004080 
00004090 
00004100 
00004110 
00004120 
00004130 
00004140 
00004150 
00004160 
00004170 
00004180 
00004190 
00004200 
00004210 
00004220 
00004230 
00004240 
00004250 
00004260 
00004 270 
00004 280 
00004 290 
00004300 
00004310 
00004320 
00004330 



p. tilt til' 3-* l ) -C5" 1 1 ( NP 3-3 ) )/tlH2 
; 4 >7 tl r ( HP3 — !5 ) 250 * 530*258 
?38 K-0 S 0 
300 DO 4'>0 1 = 1*1. 

"50 DO 351 J=3*KP3 

PlJJMIt JJ-l'-'UtJ-D-QMltJ-Z) 

351 CtJ) -Dt J)-P"CtJ"l!-0«CU“2) 
IFtllW3-l) ) 352 *400 * 35 1 

352 ir(B(NP3-l) 1353. 400*353 

353 AVH81-' AB5I t H t MPO-1 1 /B ( NP3-1 ) ) 

356 I P t AVI-13 1--SK ) 450 * 3 5 4 . 3 5 4 

35 A 0(MP3)=iH!:P3)-Q-B(NP3-2 ) 

AOO IF (B ( MP3 M '-01 *550 1 AO 1 
A01 AVI1R2=ABSF(H(I!P3) /OtMP3 1 ) 

•'03 IT . •'”"A , 'H , 32 1 550 * A50 , A 5 0 
A 30 CO A51 J---3*fiP3 

Dt J J -II t JH-R"DI J-l) 

-.3,1 itJ) =D(J) 

!FtD(i;P31 ’ 62 *500* A52 
A52 Avv:03=.' B6F t H ( IIP3 ) /D ( MP3 ) I 
'-30 IF ( SK-AVHD3 ) 500 * A53 * A53 
A3 3 CC2“-C(MP3-21 
CC3 Ci (1P3-3 1 
C( I-P3-1 ) --P- CC2-Q" CC. 

CCl=CtMP3-l) 

S=CC25CC2-ca-"CC3 
IF (S) A55 jA 5A* A55 
A5A p-r -2.0 

Q-O:: tO i-l.O) 

GO TO 436 

•A55 P--P-E ( B ( FP3 -1 ) ■'•CC2—B ( MP3 I :; CC3 ) / S 
Q=0-i t-B (HP3-1 ) : CC1 : B t MP3 ) ’'CC?. ) /S 
A56 IF(E(MP3-m 450.457*458 
A37 P.--R- 1.0 
GO TO A90 

A53 R-'ll— DtilPT ) /EtMP3~l ) 


0000A3A0 
OOOOA3SO 
OOOOA360 
OOOOA370 
00004300 
00004390 
OOOOAAOO 
0000441 0 
O000AA20 
00 00 A 430 
OOOOAAAO 
00004A50 
OOOOAAOO 
0000AA70 
00004A80 
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00004500 
0000A510 
009CA320 
00004530 
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00005550 
00094560 
OCOO 4570 
0 '■405500 
O jti 590 
54600 
04610 
i,’ 004620 
00004630 
00004640 
00004650 
00004660 
00004670 
00004600 
00005690 


no comumue 

PS=PT 

OS--QT 

PT-p 

OT =0 

IF (REV) 4° 1 ,-49? *492 

491 Sr4=SR/ie • 

492 RCV= RCV 
C-0 TO 250 

500 i F ( T ) 50 1 « 502 *502 

501 R=1.0/R 
592 MP--RP3-3 

UtMP)--R 
Vi KP I =0.0 
CGMV f NP ) =S< .. 

MP3=KP3- 1 
DO 503 J=3 » MP3 
^03 IHJ)=D(J) 

IF ( MP3-3) 300 *51.300 

550 I Ft T ) 553 s 552 » 552 

551 P=P/Q 
0=1. 0/0 

552 PP2=P/2.0 
Q!!PSQ=Q“PP2"PP2 

560 IFtC.iPSO) 5D4»5y4»553 

553 fSP=RP3-3 
U(MP)=-PP2 
O(RP-l) =-PP2 

£ -SQRTF tOJiPSOI 
V ( MPO =s 
V!MP-1)-- S 
G-0 TO 561 

354 5-50P I'F t • OliPiO) 
UP.-KP3-3 
irtP)555>556*336 
555 U(MP!--fP2-l5 
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Fig. 1-2. Double Lead Compensation of Plant 
with ,G(s) = -~r 


Fig. 13. Double. Lead Comjpfensatiori of Plant 
with G(a) = 










: '.Fig. 1-7. Double Lead' QompenWation of Plant 

with G(s)' = ^ 1 ^ 

s (s+100) 









Fig. 1-8. Double Lead .Compensation of Plant with 

with G(s) = 

<s+5)^(S+10) 


2 

TRANSIENT RESPONSE OF NONLINEAR SYSTEMS 
2.1 INTRODUCTION 

When a system has one nonlinear element that is single val- 
ued and non-frequency dependent, analysis 'of the— system is con- 
veniently accomplished using the parameter plane methods. The 
nonlinear element is represented by a describing function, which 
is a function of signal amplitude only. The describing function 
is designated as one of the parameters, a or 0. This designation 
removes the nonlinear parameter from the functions that determine 
the parameter plane curves* so that these may be plotted on the 
a-j9 plane. The M-point is located on' the a-P plane in the usual 
way, but for the case of one nonlinear element one coordinate of 
the Mrpoint is the numerical value of the describing function of 
the nonlinear parameter. For linear systems the M-point is sta- 
tionary on the cl-$ plane, but for a nonlinear system the M-point 
moves because the numerical value of the describing function is 
a function of signal amplitude. For a system with one single 
valued nonlinearity, N, where N is designated as j8, the locus 
followed by the M-point is a straight line parallel to the 0-axis. 
This locus of M-point motion can be said to start at the value of 
/3 corresponding to very small (zero) signal amplitude into the 
nonlinear element. The displacement of the M-point along this 
locus is determined by the way in which P varies as a function 

of signal amplitude, and this is determined by using the 

Constant -C and constant -a> n curves, or constant -cr and constant 
-w curves. 





the describing function of the nonlinear element. 

Previous work has shown how to predict limit cycles using 
K-point locus on the parameter plane. If this locus crosses the 
stability boundary (C=0 curve or cr=0 curve) the ifttersection of 
these curves defines the frequency of the limit cycle. If an 
amplitude -scale can be determined for the location of the M-point 
on the M-locus, then this scale is used to define the amplitude of 
the limit cycle. 

The concept of a moving M-point on the parameter plane can 
be used to calculate the transient response of nonlinear systems. 

As the M-point moves along the M-locus, each point defines both 
signal amplitude and all roots of the characteristic equation. 

This information can be used to determine the amplitude vs time 
relationship which is the transient response. Computations are 
based on Siljak's extension of some basic work by Krylov and 
Bogoliubov, and details are given in the following paragraphs. 

Assume* that the system is second order, and that the non- 
linear element is represented by its describing function. Then 
for an initial signal amplitude A q , the transient response is de- 
fined by 

X(t) = A o e CTt cos ( at + <t>) (2-1) 

where a and co are botb functions of the signal amplitude. 

*These assumptions restrict use of this method to systems in which 
a pair of complex roots dominate the transient response, and these 
systems must have low pass. filter characteristics to justify use 
of a describing function. 
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(2-2) 


CT = d(A) 

W = W(A) 

The parameter plane curves are prepared, the M-locus is super- 
imposed on them, and the describing function is used to associate 
an amplitude scale with the M-locus. Then the values of cr (A) 
and to(A) may be read from the parameter plane for any X. 

The transient response of the system from any initial dis- 
placement, A q , is determined in two steps, the first of which is 
to calculate the envelope of the transient. Assuming that 0 
(in eqn. 2-1) is zero, the envelope is defined by 

X(t) = A o " (2-3) 

which may be approximated over a short time interval by a straight 
line tangent to the exponential curve. Thus at t = 0, X=A 0 and 

from the parameter plane o-(A q ) is evaluated. Then X(t)=A Q e ° 
is approximated by a short straight line segment on the X vs t 

plane. This straight line iB terminated at t = tj and at t 1 a 
new amplitude A^ is read from the curve. Entering the M-locus on 
the parameter plane with A^ values are obtained for cr^ and to^. 

The envelope of the transient is extended from t^ to ^ with 
another straight line segment defined by X = A^e 1 , This pro- 
cedure is repeated until the envelope is defined over an accept- 
able time interval. 

As a by-product of this procedure, a) has been determined 
quantitatively as a function of amplitude and also as a function 
of time. Using the definition 
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(2-4) 


$ = J <c(a) at 
“o 

the phase can be determined at any t by graphical integration 
(i.e., evaluation of the area under the M(A) vs t curve). If 0 
in eqn. 2-1 is zero, then X(t) - 0 for §= (2n-l) (ff/2) . Values 
of t corresponding to 4> = 90°, 270°, 450°, etc., are determined by 
graphical integration, are marked on the axis of the X Vs t plane, 
•and the transient response is drawn tangent to the envelope and 
intersecting the X= 0 axis at the indicated values of t. 

The above procedures are readily applied to systems with one 
nonlinearity, and correlation with simulation results is excel- 
lent. Since such applications are elementary no illustrations 
are given here, and the study is extended to systems with two 
single' valued nonlinear elements. In general no other methods 
exist for predicting the transient response of systems with two 
nonlinear elements, so the results obtained here represent a 
significant advance in the state of the art. 

2.2 CLASSIFICATION OF SYSTEMS WITH TWO NONLINEARITIES 

When a system contains two nonlinear elements, N^ and N 2 , 
that are single valued and ’are not frequency dependent, parameter 
plane representation may be used but both a and J3 become functions 
of N^ and Nj. Computation of the parameter plane curves presents 
no difficulty, but determination of the M-locus may be difficult. 
As a result it is convenient to classify nonlinear systems accord- 
ing to the structural conditions which complicate the evaluation 
of the M-locus. The following classes are proposed: 
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CLASS 1. Identical signal excitation to both nonlinear 
elements . 

In Fig. 2-la, the signal X is the input to both nonlinear 
elements N^ and . For every value of X corresponding values 
of and N 2 are uniquely defined and are independent of fre- 
quency so evaluation of the M-locus is easy. 

CLASS 2. The input signals to the two nonlinear elements 
are related by a linear differential equation. 

In Fig. 2-lb the signal X is the input to N 2> but the 
input to N^ is X G .^(s). Thus the input to N 2 is a function of 
amplitude only, but the input to is a function of both ampli- 
tude and frequency. For a given amplitude of 'the signal X, the 
describing function for N 2 provides one unique value, but for 

each amplitude of X the describing function for N^ has an in- 

. . *. 

finite number of possible values, one for each possible value of 

frequency. As a result the evaluation of the M-locus is consid- 
erably more difficult than for Class 1. 

CLASS 3. The input signals to the two nonlinear elements 

are related by a nonlinear differential equation. 

Fig. 2-lc illustrates this class of nonlinear systems. The 
signal X is the input to N^, but the input to N 2 is X Cn^ } Cg 2 (s) ] 
where the brackets are intended to represent some functional re- 
lationship rather than a multiplication. Evaluation of the M- 
locus can be very difficult for such systems. 

2.3 EVALUATION OF THE M-LOCUS. THE DYNAMIC DESCRIBING FUNCTION. 
When a system with one single valued nonlinear element is 

represented on the parameter plane the M-locus is clearly a 
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straight line parallel to one of the coordinate axes. Thus the 
M-loeus itself is readily found hut the amplitude scale asso- 
ciated with this locus must be evaluated. For systems with two 
nonlinearities (especially Class 2 or 3) the path of the M-point 
on the parameter plane cannot be predicted by inspection. It 
can be calculated, however, using the ordinary describing function 
to define the amplitude relationships. 

To justify the choice of the describing function as a tool, 
consider the fact that parameter plane predictions of limit cycles 
are defined on the basis of a single point where the M-locus in- 
terests the stability boundary. This single point defines both 
the fundamental frequency of the oscillation and also the ampli- 
tude of this fundamental component. It is clear that the loca- 
tion of the M-point represents some sort of average value of ampli- 
tude, since the instantaneous value of amplitude varies cyclically 
during a limit cycle. The describing function of a nonlinear 
element effectively averages the response of the element to a 
sinusoidal input over one cycle of operation. Thus its use is 
clearly justified when system operation is periodic and lightly 
damped. While not so clearly justified for other operating con- 
ditions it has given surprisingly accurate results and therefore 
will be used until a better technique becomes available. 

Using the describing functions of the two nonlinearities in 
a system, a family of decribmg function curves are computed and 
plotted on the a-j3 parameter plane. When these curves are super- 
imposed on the regular parameter plane curves, the M-locus can 
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be determined. The M-locus represents the curve along which the 
M-point moves when the system is in dynamic operation, and it 
consists of the locus of all points at which -the describing func- 
tion curves and the parameter plane curves have common frequency 
intersections. We choose to call this curve the "Dynamic De- 
scribing Function Locus". The procedure and also a justification 
is as follows: 

a) Assume a constant amplitude, constant co signal at X, the 
input to one nonlinear element. Using the describing 
function compute the equivalent gain of that element; 
also compute the signal amplitude at the input to the 
second nonlinear element, and the equivalent gain of 
this second element. 

b) The two equivalent gains evaluated in (a) determine one 
point on a describing function curve on the a-P plane. 
Repetition using the same value of d> but different 
amplitudes at X determines a describing function curve 
for a constant a> signal. 

c) Repetition of a) and b) for other values of w provides 

a family of describing function curves, each curve being 
for a designated value of to. 

d) These curves are then superimposed on the usual* para- 
meter plane curves. The constant -a) describing function 

-*Curves for constant -O' and constant to are most convenient, but 
constant -C and constant W n curves can be used if it is noted that 

w = w n ./l-C 2 . 
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curves will intersect the constant -to parameter plane curves, 
and those intersections for which the (0 is the same. Define the 
Dynamic Describing Function locus. 

The nonlinear system is described by one nonlinear dif- 
ferential equation. The procedures used here effectively parti- 
tion this equation into two parts, a linear part represented by 
the parameter plane curves, and a nonlinear part represented by 
the describing function curves. Then parts are "coupled" by the 
parameters a and & which are the coordinates of both plots. If 
the system is in steady state periodic motion at a given fre- 
quency the- nonlinear differential equation of the system must be 
satisfied, so the linear and nonlinear partitions must be satis- 
fied at that frequency. This condition can exist only at the 
intersection of the common frequency curves. The points thus 
defined oh the "Dynamic Describing Function Locus" are deter- 
mined on the basis of steady state sinusoidal operation (unforced) . 
Under transient conditions the M-point moves along some locus on 
the parameter plane, and we assume that the points on The Dynamic 
Describing Function locus apply to transient operation although 
they are determined by means of steady state sinusoidal concepts. 
Experimental results indicate that this is a good assumption. 

2.4 CALCULATED AND EXPERIMENTAL RESULTS 

In order to verify the correctness and the applicability of 
the dynamic describing function’ and the graphical transient re- 
sponse calculations, specific examples of each of the three gener- 
al cases of Fig. 2-1 were investigated. The' details of some of 
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these examples, and the corresponding calculated results are pre- 
sented here. Simulation of the systems provided experimental re- 
sults which are also presented to permit comparison between theory 
and experiment. 

System 1 . _ Two nonlinear elements with identical excitation: 

The block diagram is given in Fig. 2-2. The characteristic 
equation is 

s 3 + 10s 2 + (10N X + 10N 2 )s + lOONj. = 0 (2-5) 

and it is convenient to let N^ = a, n 2 = 0. Fig. 2-3 gives the 
parameter plane plot (in o- and curves) . Since the two non- 
linear elements have identical excitation a single dynamic de- 
scribing function curve is obtained which’ is independent of fre- 
quency. However, the dynamic describing function is dependent on 
the specific numerical characteristics of the nonlinearities, and 
Fig. 2-3 contains three dynamic describing function curves (dotted] 
for three different sets of characteristics in N^ and n 2 . These 
three curves were chosen to illustrate different root variations. 

For curve 1 a real root becomes dominant early in the transient, 
for curves 2 and 3 complex roots are dominant, the system being 
moderately damped for curve 2 but going to a stable limit cycle 
for curve 3. 

Calculated and analog computer results are given on Figs. 2-4, 
5,6. It is seen from Fig. 2-4 that the dominant real root condi- 
tion cannot be handled accurately with the graphical computations. 

It is not known whether the discrepancy lies solely in the graphical 
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method which is based on complex roots, or whether the dynamic 
describing function also contributes to the errors. Research on 
this point is continuing. For the cases of Fig. 2-5 and 2-6 the 
calculated- results compare well with the computer results. 

System 2 . Tito nonlinear elements related to a common signal by 
a linear differential equation. 

The block diagram is given in Fig. 2-7, and the parameter 
plane curves With dynamic describing function curve shown dotted 
are given on Fig. 2-8. Fig. 2-9 gives the describing function 
grid needed to obtain the dynamic describing function curve. To 
obtain the grid of Fig. 2-9 the point A q on Fig. 2-7 was chosen 
as a reference point, and at each value of to the amplitude of the 
(assumed) sinusoidal signal at A q was varied to obtain the vs 
Nj values for a constant to curve on Fig. 2-9. The dynamic des- 
cribing function curve on Fig. 2-6 is obtained by superimposing 
the parameter plane curves of Fig. 2-8 on the describing function 
net of Fig. 2-9 and locating intersections of constant to curves 
of the same to value. 

Limit cycle predictions of the dynamic describing function 
curve on the parameter plane agree with analog computer simulation 
results. In addition Figs. 2-10, 11,12 compare predicted transient 
response with simulation results. 

Additional checks were run using different values for the 
deadzone and saturation limits m the two nonlinearities, ’ but 
the detailed data is not given here. In general the predicted 
and simulated results were in good agreement except when a real 
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root became dominant during the transient response, in which case 
the frequency of the oscillatory component was usually predicted 
with reasonable accuracy, but amplitudes were not, nor was the 
total response time due to the influence of this real root. 

The calculations and simulations were also repeated with the 
nonlinearities interchanged (i.e., in Fig. 2-7, becomes a 
saturated element and N 2 a dead zone element) . Using the same 
techniques the results obtained were always in agreement with 
about the same degree of accuracy and with the shortcomings as 
previously noted. 

System 3 . Two nonlinear elements related to a common signal by 
nonlinear differential equation. 

The classification described as System 3 can contain a wide 
variety of combinations of linear and nonlinear elements, of 
which the parameter plane method may be applicable to only a 
small subset. A specific system which belongs in this class is 
shown in Fig. 2-13. The characteristic equation of this system 
is 

s 3 + 3s 7 + 2s + 40KN^ (H^ + gt^) (2-6) 

where Nj = N a + for the hysteretio nonlinearity, and we define 
a = N^N a ; (9 = N^H^. The parameter plane equations are still appli- 
cable and the parameter plane curves can be computed. For the 
purposes of this study only C = 0 curve was calculated, and only 
the limit cycle predictions were checked. The describing function 
net is required, and in this case relates the NjN a and pairs 

to the common signal at A on Fig. 2-13. The results of these 
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computations are given on Fig. 2-14, which shows the E = 0 curve 
from the parameter plane equations and the describing function 
net for the case where K = 0.15. Only one point is defined on 
the dynamic describing function curve, and this is marked on 
the E = 0 curve at the point where the a) value on the £ = 0 curve 
is the same as the value of the constant w describing function curve 
passing through that point. This defines the frequency and ampli- 
tude of the limit cycle, and the results agree with simulation 
results. 

Note that a change in the value of K changes the differential 
equation of the system, thus requiring a new set of curves. Re- 
sults were obtained with other values of K and again the predic- 
tions agreed with simulation results. 

2.5 COMMENTS 

The results obtained thus far indicate that the parameter 
plane is a useful tool m predicting the stability and response 
of nonlinear systems. The accuracy available is only fair, but 
is more than adequate for many engineering applications. The 
transient response predictions - in particular for systems con- 
taining two nonlinearities, - are better than are available with 
any other method. 

The graphical presentation of the dynamic describing function 
curve on tne' parameter plane is potentially a valuable design 
tool, it indicates at a glance the range of variations of the 
roots, and thus permits prediction of a desired location of the 

describing function curve, which in turn implicitly defines the 
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required characteristics of the nonlinear element. Further 
research is required in this area. 

The technique becomes inaccurate when the transient re- 
sponse is influenced by more than two complex roots. Again more 
research is required to evaluate this situation. 

It is too early to assess the true value of studying non- 
linear systems on the parameter plane. Without question it does 
make possible many types of analyses that are not readily avail- 
able otherwise. However, the limitations of the technique are 
not clearly defined, and it obviously is important to know under 
what conditions the methods are not applicable, or should be 
applied with care. 
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Fig. 2-1. General Classification of Control Systems 
with Two Nonlinear Elements. 














Fig. 2-3. Dynamic Describing Function Curve on 
Sigma-Omega curves. 
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pig. 2-7. Block Diagram of Fourth Order System with 
Two Nonlinear Elements. 











Pig. .2-8. Dynamic Describing Function Curve on 
Sigma-Omega Curves. 
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Fig.. 2-12a. 


Calculated Transient Response for A - 3 volts." 



Fig. 2 -12b. Analog Computer Transient Response" for A' =■ 3 Volts. 
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Fig. 2-13. '.Block Diagram .of Test Example. 
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3.1 Introduction. 


CHAPTER III 

ASYMMETRICAL NONLINEAR OSCILLATIONS 


In certain classes of nonlinear systems, oscillations may 
consist of a limit cycle superimposed on a constant or slow-varying 
signal. These oscillations are referred to as asymmetrical oscil - 
lations since the center of the limit cycle is shifted according 
to the corresponding value of the constant or slow-varying signal. 
In general, asymmetrical oscillations may occur when the input- 
output characteristic of the nonlinearity in the system is not 
symmetrical about the origin, or when the system is subject to 
forcing signals. When the nonlinear characteristic is asymmetric, 
the output of the nonlinearity may contain a constant term even 
though the corresponding input is a single sinusoidal wave. If 
the nonlinear characteristic is symmetric, asymmetrical oscil- 
lations can arise whenever the system is subject to forcing input 
signals. Evidently these oscillations may take place at certain 
points of the system if both conditions are present. Before the 
analysis of asymmetrical oscillations in the parameter plane is 
presented, the previous work and results in’ considering these 
oscillations and related problems are reviewed. 

It has been shown first by MacColl C 3 . 1 D that the introduc- 
tion of an external sinusoidal signal at the input to an on-off 
servomechanism yields a system that behaves like a linear one for 
small inputs superimposed on the sinusoidal signal. This pheno- 
mena has been later investigated under various names, such as 
"dither effect", "signal stabilization", etc. Asymmetrical non- 
linear oscillations has been found by a majority of authors as 
the most appropriate term for the mentioned phenomena. 
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In analyzing a carrier-controlled relay servo, Lozier [3.2] 
has introduced an idea to accomplish the linearization of the re- 
lay by a limit cycle existing in the system and without an ex- 
ternal signal. This idea has been further developed by several 
authors [3. 3-3. 9] and a detailed treatment of the problem has been 
given by Popov and Palitov [3.8], On the other hand, the external 
signal application has been developed by Loeb [3.9] and Olderiburger 
with his associates [3.10-3.12]. The latter introduced the name 
"signal stabilization" to indicate that the nonlinear system is 
stabilized in the state of sustained oscillations with suffici- 
ently high frequency. The stabilization is actually a consequence 
of the linearizing effect discovered by MacColl. The concept of 
signal stabilization has been extended by Sridhar [3.13-3.14] to 
the case of a nonlinear system which has one single-valued non- 
linearity m the loop, and the stabilizing signal is a stationary 
random process with a Gaussian distribution and obeys the ergodic 
hypothesis. 

The above defined problem can be treated by dual -input de- 
scribing functions as proposed by West [3.15]. This approach has 
been significantly simplified by Boyer [3.16] as outlined by 
Gibson [3.17]. A similar approach is used by Gelb and Van der Velde 
[3.181, and significant results have been obtained by Atherton and 
others [3.19-3.20] who made a comparison of the utilized concept 
with the Tsypkin method [3.2l], 

The study of asymmetrical nonlinear oscillations has been 
extensively performed in tho analysis and design of a large class 
of plant adaptive control systems. This class of system is 
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sometimes called the limit eyeliner adaptive systems because of 
the fact that the existing limit cycle is used as an identifi- 
cation signal. Some of the references on this subject are listed 
here [3.22-3.26]. A majority of the authors proposed an external 
sinusoidal signal for identification. More recently, Gelb and 
Van der Velde [3.18] have examined .to a limited extent and in 
a quantitative manner the properties of self-oscillating adaptive 
systems which have several advantages over the external adapta- 
tion, such as simplicity, cost, reliability, etc. The following 
analysis of asymmetrical nonlinear oscillations in the parameter 
plane can be applied directly to self-oscillating adaptive sys- 
tems. 

In the following developments, the asymmetrical nonlinear 
oscillations are analyzed in the parameter plane [3.273. The 
control systems with asymmetrical nonlinear characteristics are 
considered to determine stability and sustained Oscillations. 

The same type of oscillations is investigated m nonlinear con- 
trol systems subject to constant reference and perturbing input 
signals. The procedure is further extended to the analysis of 
systems with slow-varying input signals. In 'this case, it is 
shown how a nonlinear characteristic can be- mo'dified for the 
slow-varying signals. The presented analysis is performed with 
respect to both input signals and the values of adjustable sys- 
tem parameters. The analysis procedure is illustrated by examples 
in which multiloop feedback structures with several adjustable 
parameters are considered. In addition, various nonlinear 
characteristics are used in either the forward or the feedback 



3-4 


path. The obtained results are checked by. computer simulations 
which indicate a sufficient accuracy of the presented procedure. 
3.2 Basic Developments 

Consider a nonlinear system described by the nonlinear dif- 
ferential equation 

B(s)x + C(s) F{x,sx) = H(s)f, s ^ JL (3.1) 

where B(s), C(s), and H(s) are polynomials in s and the degree of 
the polynomial B(s) is greater than the degree of the polynomials 
C (s) and H(s) . The function F(x, sx) describes the nonlinearity. 
Function f = f(t) is a forcing signal, which may be either a 
reference input or a perturbing signal, and it is assumed to be 
a constant or a slowly-varying function of time. 

As a first approximation, the steady-state solution x = x(t) 
of equation 3.1 which represents the input to the nonlinearity, 
is assumed to be 

x = x° + x* (3.2) 

where x° = x°(t) is either a slowly-varying function of time or 
is constant, and x*, which is' 

x* = A sm 0, 0 = C<t + 0, (3.3) 

represents the periodic component of the solution x(t). Since 
6 ip (3.3) merely corresponds to a shift m t, one can put 0=0 
and use x* = A sin Clt. 

The forcing function f(t) is considered as a slowly-varying 
function of time if it can be assumed approximately as constant 
ova any cycle of the periodic component x*; i.e., 

|f(t+T) - f (t) | < - |f(t) 1 


(3.4) 
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where the period T = 2ff/CI. In the frequency domain, equation 3.4 
means that the frequency SI of the periodic component x is much 
greater (practically ten times or more) than the highest frequency 
of the slowly-varying component x°. In this case, no harmonic 
relation between the components x° and x* nonlinear system sub- 
ject to forcing signals, such as jump-resonance, generation of 
subharmonics, etc., cannot take place. The forced nonlinear os- 
cillations for which the condition (3.4) is not satisfied neces- 
sarily, are considered in other works. 

Under the condition (3.4), the values of x°. A, and 0, which . 
appear in the solution x = x° + A sin fit, are slowly-varying 
quantities in time. This enables the extension of the conven- 
tional harmonic linearization in which the describing- function 
is defined for the signal x = x° + x* as an input to the non- 
linear element. Thus, the nonlinear function F(x,sx) is approxi- 
mately expressed by the first terms of the Fourier series as 

N 

F(x,sx) = F° + N-jX* + -j^sx* (3.5) 

where 

.2ir 

F° = j F(x° + A sin 0, AO cos 0) d0 (3.6a) 

o 

,, P 2ff 

W 1 = 7TA J + a sin 0, a CJ cos 0)sin 0 d0 (3.6b) 

i r 2,r 

Nj “ J F(x + A sin 0, A Cl cos 0)cos 0 d0 (3.6c) 

and 0 = Cl 1 1 

As can be seen from equations 3.5 and 3.6a, the component 
F° of the output of the nonlinearity F(x,sx) is not considered 
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zero as was the case m the analysis of symmetrical nonlinear 
oscillations presented in the previous chapter. This results 
from the fact that either the nonlinear function F(x,sx) is not 
symmetric or the system is subject to an external input signal, 
or that both facts are present in the system. 

According to equations 3.6, all coefficients F°, EJ^, and N 2 
are generally functions of x°, A, and 0, i.e., 

F O +F°(x°,A,0), ^ , K 2 = W 2 (x°,A,fi) (3.7) 

For a majority of the nonlinear functions F(x,sx) encountered in 
practical applications, the above functions (3.7) are obtained 
once and for all. 

By applying the linearization of the function F(x,sx) given 
in equation 3.5, the solution x = x° + x* of (3.1) can be ob- 
tained by considering the following linearized differential 
equation 

B(s) (x°+x*) + C(s) (F°+NjX* + -$■ sx*) = H(s)£ (3.8) 

instead of equation 3.1. If x°. A, and 0 are slowly-varying 
functions of time as a consequence of the same property associated 
with the forcing function f, equation 3.8 can be rewritten as two 
simultaneous equations corresponding to the slowly-varying sig- 
nal x° and the periodic signal x* as follows: 

B(s)x° + C(s)F° = H(s) f (3.9a) 

. N 

B(s)x* + C(S) (H.X* + -ff- SX*) = 0 (3.9b) 

Equations 3.9, however, cannot be solved independently since they 
are related to each other by the nonlinear equations 3.7. This 
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fact indicates that the applied linearization preserves the 
essential feature of nonlinear systems and that the superposi- 
tion principle from linear analysis is not valid. 

An analytical solution of equations 3.9 is difficult to 
obtain since F° in(3.9a) is usually a trancendental function 
with respect to x°. A graphical procedure is presented for 
solving equations 3.9 in the parameter plane. A necessary con- 
dition for equation 3.1 to have a solution x(t) close to 3.2 is 
that the characteristic equation 

KL 

. B(s) + C (s) (N x + ~s) = 0. (3.10) 

corresponding to the linearized differential equation 3.9b, have 
a pure imaginary root s = j£i. 

By using the parameter plane approach, equation 3.10 can be 
solved for a and 0 as 

a = a(i^ (3.11) 

0 = P(O) 

where a and 0 are and H 2 or some other system adjustable 
parameter. Equations 3.11 represent the £ = 0 (or C = 0) curve 
for which s = jO. The £ = 0 curve determines the stable region 
in the a0 plane in the usual manner. After the stable region is 
found, the loci of points M(a, 0) "are plotted according to the 
variations of a and/or 0 representing and/or N 2 . — The M loci 
incorporates the additional variable x°, and a family of the 
loci should be constructed for different values of x°. Then the 
stability of the nonlinear system is determined by the relative 
location of the £ curve and the M loci and the limit cycles are 
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found at their intersections.. The stability of r the limit cycles, 

is determined in the usual manner. This part of the solution 

pfocess will be best described by the examples that follow. 

The presence of a limit cycle in the system can modify the 
nonlinear characteristic for the slowly -varying input signal. ■ 

In order to determine the modified characteristic, the intersec- 
tions of the S = o curve and the M loci are considered' to evalu- 
ate the amplitude A and the frequency Q of the limit cycle as' 
functions of the slowly-varying component x°; i.e., 

A = A(x°), a = fi(x°) (3. ,12) 

These functions, when substituted into the function F° (x°,A, fi) 
yield the modified nonlihear characteristic for the slowly- 
varying signal 

>° = >Mx 0 ) (3.13) ' 

The function iji (x°) is continuous' in a limited range of x°‘, which 
indicates the smoothing effect due to the presence of the limit 
cycle. 

Substitution of equation 3.13 into equation '3.9a gives 

B(s)x° + C(s) =H(s)f • ' (3.14). 

Equation 3(14 is a nonlinear differential equation in' x°, which 
can be solved graphically for x° after the function <! >(x°) is 
obtained. This, in turn, yields the related values of the 
functions A(k°) and fi(x°) of equations 3.12,' and 'the solution 
x . = x ° +, A sin fit is thereby determined. 

The function $ (x°) is a continuous function of x° and it can 


be assumed approximately linear for small variations of x°. Then 
the stability problem related to equation 3.14 can be solved by 
known linear methods. If it is regarded as a nonlinear function 
of x°, it can be linearized by harmonic linearization and the 
results of the previous chapter can be applied. 

It should be noted here that the same parameter plane pro- 
cedure can be used when the right side of equation 3.1 has more 

than one forcing, function; i.e., the right-hand side is expressed 
r • * 

by E H. (s)f. . The solution x, however, must be found by con- 
i=l 11 . ’ 

sidering all existing inputs simultaneously since the super- 

. * T 

position principle of linear analysis is not valid. Further- 
more, if the . polynomial H(s) of equation 3-1 can be factored. in 
the form sH^(s), the procedure applied to the case in which the 
rate sf of the .function f is considered as a slowly- varying sig- 
nal; i.e.,'' |sf(t+T).- sf (t) J . 

The presented graphical procedure can be extended to. non- 
linear 'Control sys'tems with two .nonlinear functions F^(s) and 
Ij'M » whereby the following nonlinear differential equation is 
investigated; 

B(s)x+. C(s) f^(x) + D(s) F 2 (x) = H(s)f. (3.15) 

' In this case, a procedure similar to that given in Section 
can be extended to determine' the solution x = x° + x*. 

The general procedure outlined in this section is modified 
depending on the actual problem involved. These problems may be 
divided into three ma^or groups; asymmetrical nonlinearities; 
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constant forcing signals; and slow-varying signals. In the 
following, each group is considered separately. 



In equation 3.16a, which corresponds to equation 3.9a, there is 
no forcing slowly-varying function (f = 0) , and in the steady- 
state solution x = x° + x*, the x° is constant and hence s is 
replaced by zero in B(s) and C(s). 

In practical situations, B(o) or C(o) can be zero. Also, 
the nonlinearity, m. the system is often .described by a single- 
valued function F(x) and N 2 =0. 'Thus, an adjustable parameter 
appearing in B{s) or C(s) can be chosen as one of the axes in 
the parameter otfi plane, while the other axes is related to the 
describing function coefficient N^. Some of these situations are 
discussed in the following examples. 

Consider a feedback control system with the block diagram 

.of Fig. 3.1 m which the transfer functions are 

K K 

G 1< S )= K 1- G 2< e ) = sTs+lj ' G S = s+2' G -1^ s )“ K - 1 s ' ( 3 * 17 ) 
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She nonlinearity n has the form shown in the 'Upper left corner 
of Fig. 3.2. 

Equations 3.16,' for the system under investigation, have 
,the form 

F° = 0 (3118a) 

{s (s4l) (s+2) + [10,1^3(342) 4 = 0 (3.18b) 

where, according to the function F(x) of Fig. ,3.2 and -.equations 


3.6, . one has 



ilrmlc Il+ffilc a . sin 3^. 
,2 V .A 



(3.X9a} 

(3.19b) 


(3.19c) 


and x = x(t) is lie input signal to the nonlinearity n as indicated 
in Fig. 3.1. . ; 

The characteristic equation of equation 3.18b is 


s (s+1) (s+2) 4 [k 2 K_ 1 s(s 42) 4 KjiyKglNj = -0 (3.20) 

By denoting K2 K -l K l = a and = P , ' the £ = 0 curve is . 

obtained, as ' 

a = | - 2) ' (3.21) 

J3 •= |^(S7 2 4 4) 

and the stable region is determined m the a (3 plane in the usual 
fashion as shown in Fig. 3.2. 

From equations 3.18a and 3.19a, one obtains' 

x 0 = A cos — — (3.22), 

14m ' , ' 
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and of equation 6.19b becomes 



By using equation 3.23 and the expressions a = k 2 K_jN-^., 

P - thre,e M- loci (a), (b) , and (c) , are drawn in Fig. 

3.'2. They correspond to the parameter values m = 0.5, c = 1, 

K 2 = 1 and (a) K^, K_ x = 0.125? (b) ’K-jK., =' 8.39, K_ 1 = 0.28? 

(c) = 26, K- ^ = 1.75. -The stable asymmetrical oscillations 

are found at the poirit ^ and m 2 'where the M loci (a) and (b) • 
intersect the -£ = 0 curve. The amplitudes of the. oscillations are 
approximately A^ = O’. 85 and A 2 = O'. 8, which is, read from the 
M loci (a) and. (b) at the--intersections M-j, and m 2 - The corres- 
ponding frequencies 0^ ='1.5 'and f? 2 = 1.6 are indicated on the ' 
f = 0 curve. The related’ values of x° in the solution x = ^ > 4CX'sin !3t 
is calculated for each point and M 2 using equation 3.22, namely, 
x° = -0.42 and x 2 = -0.39. 

Xn Fig. 3.3, the solution x^ = 0.42 4 0. 85 'sin '1. 5t for the • 
.case (-a) is shown as obtained by a digital computer 'simulation. 

The calculated results are sufficiently close to that obtained, 
by the simulation.' -From Fig. 3.3, it can be seen that an initial 
'condition x-^to) = 4.25 is used and the variable Xj (t) approached 
a stable limit cycle. That the limit cycle is stable and will be 
reached by x^(t) starting from x^(o) =4.25 can be concluded from 
the relative location of the £ = 0 curve and the M locus (a) , as 
explained in the preceding chapter on the symmetrical oscilla- 
tions. The additional component x° of the solution x(t) does 
not alter the stability analysis of .the oscillations. 
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An analog computer simulation of the case (b) gives the ' 
solution x 2 = -0.39 + 0.8 sin 1.6t as shown in Pig. 3.4. A suf- 
ficient accuracy is indicated. .The initial condition x 2 (o) = 0 
and x 2 (t) reached a limit cycle. This could be concluded from 
Fig. 3.2 as previously noted. 

It is of particular interest to consider the case (c) of 
Fig. 3.2. The M locus (c) is tangent to the = 0 curve and cor-, 
responds to the ratio a/p ~ KjKj/K^ =, 14.8. If this ratio is 
higher than 14.8, then there is a limit cycle as shown by cases 
(a) and (b) . On the other hand, if this ratio is less than 14.8, 
the entire M locus is situated in the 'stable region and 'the cor- 
responding system is always stable. The tangent case (c) : 
m = 0.5, c = 1, K 2 = 1, Kj_K 3 = 26, K_ 1 = 1.75, is simulated on 
a digital computer and the obtained solution x 3 (t) is shown in 
Fig. 3.5, which indicates that the system is stable. 



When the forcing signal at certain points of a nonlinear 
system is constant, the solution x = x° + A sin Ot(if it exists) 
will have x°. A, and £2 as constant values. To determine these 
values, note that the equations to solve in the presence of a 
constant forcing signal f have the form 

B(o)x° + C{o)F° = H(o)f° (3.24a) 

N, * 

[b(s) + .C(s)H 3 + s) 3x =0 (3.24b) 

In general B(o), C(o), and H(o) are constants different 
from zero, and the solution, procedure is somewhat more compli— 




Syfcl a block difudvna 
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(3.28) 
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The variation of the M point due to the function Nj • Mj (x°. A) 
given «■ o 

Nj - k - f (arc sin + are sir. °^ X 

4 A - (2^)* 4 2*!L. A - (^J 5 ^)*. A >0 ^ |x°l 

(3.29) 

la plotted in Pig. .7c. (The expression (3.29), corresponds to 
the nonlinearity c Pig. 3.7a). In order to find a solution 
x » x° ♦ x* of equations 3.27, the parameter k is assumed equal 
to one, and the function P°(x°, A) is plotted in Pig. 3.7d by 

r ° . Ji - (S^V - A - * kx° -t 

• iCDIarc sin - are .In ^J 2 -) - 

- x°(arc sin + arc sin ^-)3* ***> ♦ lx°| 

(3.30) 

For T_ t - 0.04. the ooint Mj(0.04r 14.3) correaponda to a 
solution x - x° ♦ x* which will have 0-12 rad/aec as indicated 
on the curve C - 0. If - 2' . from Mj it follows that 
N 1 - J/K a ■ 0.715. This valuo of Nj detonninea the relation- 
ahip betvoon the values of x° and A for a possible solution x. 

This relationship, expressed as a function A - A(x ). can be 
graphically obtained from the diagram Nj - Njtx 0 , A) by plotting 
the straight lino PjP 2 corresponding to the value Nj - 0.715. 

The function A • A(x°) represents the solution of equation 
3.27b only. Tho pair of valuos (x°. A) which enter into the 
actual solution of equation 3.27, in rcplotted on the diagram 
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r° - r°(x° # A) of Fig. 3.7d into tho curve P{P£. Suppose t* t 
the constant perturbing signal has a value of f° • 11.75r then 
aquation 3.27a determines tha straight lina f° • 11.75 plotted 
in tha diagram r° • P°(x°. A) of Fig. 3.7d. Tha intersection R 
of that straight line and the curve PjP£ gives the pair (x°. A) 
of the solution x(t) which satisfies equation 3.27 simultaneo :<jly. 
At this point R, the values are x°/D « 1.35 and A/D ■ 1. Tho 
same values are obtained at tho point 0 on the diagram Bj ■ H^(x°,A) 
and the solution x • x° ♦ A sin Ot of equations 3.27 is found. 

If o - 1, it is x • 1.35 ♦ sin 12t. Mote that the same solution 
is obtained if the point of Fig. 3.7b is considered save that 
the frequency 0 is lower (approximately 0- 5.5 rad/sec) . 

Simpler situations may occur if one of the values B(o) or 
C(o) is zero. To illustrate, considor the nonlinear system of 
Fig. 3.8. The transfer functions are 

c i<s> - *r C 2 U) " rrrnv C J a) " alT* G -1 U) " K -l* 

(3.31) 

and tho nonlinearity n in the system is given by the function 
F(x) of Fig. 3.9. The Input to tho system is the reference con- 
stant input signal r(t) • r°. 

The nonlinear differential equation describing the above 

system is 

rs(s4l){s*2)»K 1 K 2 K 3 )x*K 2 K_ 1 s(sf2)F(x) - K 1 K 2 (s+2)r° 

(3.32) 


which may he rewritten according to equations 1.24 as 
K 1 K 2 K 3 x ° " 2t ° 


(3.33a) 





Tin. 3.11 • Function *^(A,X°) 


\ 


LZ-Z 
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[■(•+1) (•♦2) ♦ KjK 2 K 3 ♦ K 2 K_ 1 «(*+2)N l ]x*- 0 (3.33b) 

Tho characteristic equation of the equation 3.33b la evidently 

(3.34) 


e(s+l) (a+2) * ♦ K 2 K_ 1 «(b42)N 1 - 0. 

By denoting 


(3.35) 


the parameter plane diagram ia plotted in Pig. 3.10 in the usual 
fashion, tho function Nj • M 1 (A,x°), which appears as a varia- 
tion of a in the point M(ar S) is plotted in Pig. 3.11 by using 
general formula 3.6b. 

Prom equation 3.33a, one can derive tho following relationship 
between the input r°, the constant term x°, and the parameter 
0 • ky ^ 

Sf - 11- (3.36) 

x°/S 

where S is the parameter of the nonlinearity P(x) of Pig. 3.9. The 
function Sfi given in (3.36) is plotted in Pig. 3.10. 

Now, by using Pig. 3.10 and 3.11, it is possible to determine 
the sustained oscillations and their stability for various values 
of system parameters Kj, K 2 , K^, K_ 1# S. k, and the input r°. For 
example, if Kj - 1, Kj - 10, Kj • 1.75, - 1, S - 1, k - 1, and 

r° - 1.1, then tho solution of equation 3.33 is determined by the 
values x° - 1.2, A - 0.3, and 0 * 2.1 rad/sec to be approximately 


x - 1.2 ♦ 0.3 Bin 2. It (3.37) 
Por a given value of P • Kj ■ 1.75, r° • 1.1. and 8*1, the value 
of x° • 1.2 i* read from the left part of Pig. 3.10. Then the 
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value of KjK 2 0 - 17.5 determines the point H(1.2> 17.5) on the 
C • 0 curve vhero 0 - 2.1 rad/soc. At this point. K,K_jO " 
which giv*B N x - o - 0.12. Pig. 3.11 is csed to evaluate the 
amplitude A • 0.3 from the curve x°/S • 1.2. The value A - 0.3 
is read directly from the diagram Nj (A, x°) of Pig. 3.11# since 
K ■ $ ■ 1 are the parameters of the given nonlinearity in Pig. 

3.9. 

The solution (7.37) i- stable since an increase in the 
amplitude A causes the point M to move into the stable rogionr 
while a docrease in the amplitude A places the point M inside 
the unstable region of the parameter piano (Fig. 3.10). It is 
of interest to note that if the produce where 0 • is 

such that it is less than 6.4, the system is always stable since 
there isno intersections of the variation of the M point and the 
C ■ 0 curve. 

The above solution (3.37) is checked by computer simulation 
to obtain the curve on Pig. 3.12. The accuracy of the calculated 
solution is sufficiently high and. for calculated values of x°, A. 
and C, is approximately 10*. On the other hand, the computer 
solution indicates a distortion of the assumed solution 
x - x° ♦ A sin Ck Which is due to the higher harmonics present 
in the actual solution. 

3.5 5 lowly-varying Signals 

in this section, the problem of linearizing a nonlinear 
system by a high-frequency lladt cycle Is considered in more 
detail. The objective is to determine the conditions under which 
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■uch a linearization ia possible and than to comtnict the 
linearized characteristic of the nonlinearity. Thia lineari- 
zation has several practical aspects discussed in Section 3.1, 
which are based upon a general property of the linearized synteai 
that, for a United magnitude of the reference sigr.-l, behaves 
like a linear system. Therefore, results of the nonlinearitioe, 
such as dead-zone, hysteresis, backlash, etc., are eliminated. 

The procedure to achieve this will be best illustrated in the 
following examples. 

Consider the system on Pig. 3.13 with the transfer functions 

V> - K. 0j(.) - o 3 U>- - K_ t 

(3.38) 

and the nonlinearity n as shown in Pig. 3.14. The input to tho 
system is a slowly-varying reference signal r - r(t). 

The equation which describes the system is 

Cs(s+1) (s 2 *0.eu-l)+K 2 K_ 1 .A(a+l)]x ♦ KjKjK^Pix) - KjKjO(»*l)r 

(3.39) 


where the signal x - x(t) is the input to the nonlinearity. Equa- 
tion 3.39 can be rewritten in terms of equations 3.9 as 

t.(iwl) W’-flCjlCjXjP” - KjKj.ls.llr 

UKjKjKjH, X* - 0 (3.40) 

The characteristic equation of the necood equation 3.40 ia 

is(li.l) (« J *0.B».I).K,K_ ! M»«U • KjKjKjMj - ° (J.41) 

Substituting KjK.j " **• *1*2 K 3 B 1 “ **>d ■ " 3^ into equation 

3.41, ono obtains the parameter plane equations of tho C - 0 curve 
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“ - 1 ' 6 - 1 (3.42) 

fl - o.e 0(0 + i). 

The curve ( ■ 0 ii plotted in Fig. 3.14. The variation* 
of the M point aro plotted also in Fig. 3.14 according to 

- if A - <f> 2 - * ’ i- 0 ' < 3 -"’ 

Tho system parameters Kj » 1, Kj • 12.5, Kj • 10, K_j • 1 
result in the point M 1 (12.5; 45). If c • 1, this point M x gives 
> J/KjKjKj • 0.36, and the straight line P^Pj is plottod on 
tho diagram of function “ Nj(x°, A). After the diagram 
p° . r°(x°. A) is plotted in Fig. 3.15 using 

F° . ^ arc sin ^ . A » |x°l < 3 -“> 

the replotting of tho straight line PjP 2 on ti ' e <,la ^ rBIB 
yields the function (x°) of Fig. 3.15. Tho replotting procedure 
is the same as that used in the previous section: i.®.. for each 
pair of values (x°,A) read on tho straight line PjPj' the cor- 
responding pair exists in tho diagram F°(x°, A), Which determines 
one point on the curve <t>(x°). 

Function e(x°) of Pig. 3.15 is smooth and represents the non- 
linearity for the slowly-varying signal x°. Tho shape of &(x ) 
explaina tho smoothing effect of the high frequency limit cycle 
which has a slowly-varying amplitude, the valuo of Which is 
located between tho points Qj and 0 2 on the A axis of Fig. 3.14. 
Tho frequency 0 is approximately constant and has tho valuo 
P- :.7 rnd/nec. According to 4(x°), the smoothing effect of tho 
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limit cycle is present under the condition that !x°T 2.25. Por 
small values of x°, it is possible to consider t(x°) • Kx° whore 
K - const. Thon tho stability of the system with respect to 
s lowly- varying signals may be investigated by well-known linear 
methods outlined in Chapter II. In the specific example, tho 
equation of interest is 

.(•♦1) (. 2 40.8»»1) ♦ KjK ,.(..1) ♦ K KjKjKj - 0 (3.45) 

Finally, it is to be noted that for the smoothing effect to take 
place, the amplitude A should be A * |x°|, as stated in equations 
3.43 and 3.44. 

The results of the above analysis aro checked by simulating 
the system on an analog computer. Three cases are considered. 

In Pig. 3.16, the input to the nonlinearity x ■ x° ♦ A sin ft: 
and the system output x » x(t) are shown when the input signal 
is r » sin O.lt. The obtained computer solution agrees with the 
prediction. The output c(t) exhibits a smaller amplitude limit 
cycle with the same frequency. Whon tho input amplitude is in- 
crease.* five times, the diagram of Pig. 3.17 is obtained. This 
change increased x°, but the amplitude A remained almost the 
same. The frequency Cl did not change. Similar renultM occurred 
when the input amplitude increased ton times except that the 
amplitude A became slightly smaller, which agrees with the dia- 
gram of Pig. 3.14. The third case is given in Pig. 3.16. It 
should be noted from these computer solutions that the output 
sigiial c(t) represents the input signal r(t) except for the 
superiasHisod limit cycle. It can be oliminatod by introducing 
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sufficient filtering in the block G y lm) of the system of Pig. 3.13, 
or by readjusting the system parameters to obtain a higher fre- 
quency limit cycle. 

If the values of the system parameters are chosen so that the 
operating point is m 2 (21. 2x 120) of Pig. 3.14, the frequoncy of 
the limit cycle becoties higher. However, tho corresponding range 
of variations of x° is decreased to |x°| < 0.7, together with tho 
range of the amplitude A which is botvoon 0^ and Q^. This indi- 
cates that the presented procoduro is convenient to apply when 
the system parameters an' operating conditions ore changed. 

If tho nonlinearity n is changed in the system of Pig. 3.13 
by introducing a considerable dead zone D, a diagram of Pig. 3.19 
is obtainod. Tho variation of the M point is calculated by using 
equation 3.6b for the given nonlinearity of Pig. 3.19. Two cases 
should be considered separately} i.e., 

Nj - A >|x°l*D (3.46.) 


"l ' if ^ 




U°|-D « A c|x°|+D 


(3.46b) 


and the diagram Nj(x 0 , A) is shown in Pig. 3.19. By using equa- 
tion 3.6a, tho corresponding diagram P(x°, A) of Pig. 3.2C is 
plotted according to 

P° ■ £ (arc sin ♦ arc sin *"j^) • A >lx°| + D (3.47a) 

P ° • § ♦ nrc nin * P )»in x°. |x°l-D < A v |x°Ud 

(3.47b) 


» \ 
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If the points and h 2 arc '.tosan in Fig. 3.19 as operating 
points, the roplotting of the straight llnea PjPj and P ? P 4 results 
in the tvo linearized characteristics a and b of Pig. 3.20. re- 
spectively. They are cot. tructed for the values of nonlinear 
parametera c - D - 1. As can be soen from Pig. 3.20 the dead 
zone is eliminated as far as the slowly-varying signals are 
concerned. Por this to take place, it is necessary to choose 
operating conditions such that equation 3.47a is valid. This 
means that the amplitude A of the lis>lt cycle must be greater 
than |x°|*D. Otherwise the linearized characteristic t(x°) does 
not go to zero when x° » 0 since P° does not go to zero for 
x° - 0. This is indicated in Pig. 3.20 whereby P° • 0 for x° • 0 
and the dead zone is eliminated. 

By the outlined technique, it is possible to eliminate the 
hysteresis and backlash in systems with multi-valued nonlineari- 
ties. The linearization yields a single-valued function e(x°) 
which is linear in a certain limited range of values of the vari- 
able x° about the origin. To illustrate this, consider a non- 
linear system with the block diagram of Pig. 3.21 and the trans- 
fer functions 

V > • '- i '* 1 * *-»■ < 3 - 46 ' 

The nonlinear function P(x) of the nonlinearity n l« given in 
Pig. 3.22. 

The equation describing the system is 

s(s4l) (s+2)s ♦ (K^jB^JKjPfx) • 0 
Aitox harmonic linearization of 3.49, the corroeponding 


0 . 49 ) 
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characteristic equation Is 

N 2 

s(a*l) (s+2) ♦ K 1 (K_ 1 s+l) (Nj ♦ -prs) - 0 
If K. - 50, K . - 1 


(3.51) 


and s ■ jO, one obtains the C - 0 curvo as 

-A* 

9 - in. (3.52) 

The curve is plotted in Pig. 3.22. On the same plot, the varia- 
tion of the point M(Nj; N 2 ) is constructed according to 

N i ■ fx (• / i-( B ir-) J * A-(fi~^)*) a : d + U°l 


■s 


(3.53) 


and tho nonlinearity F(x) of Pig. 3.22 for which c - D - 1. Prom 
the intersections of the C • 0 curve and tho variation of the M 
point, one can determine the amplitudo A and the frequency 0 
as function of x°j i.e., 

A • A(x°) 


<x v > 


Then, by using the expression 
D»x° 


(3.54) 


p“ - *(«rc Bin - »rn .in ^-1. A > D ♦ lx°| (3.55) 

for c • D ■ 1 , a family of curves with constant amplitudo A is 
plotted on Pig. 3.23. If the first equation 3.54 is mapped onto 
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Fir. 3.2* - Computer solution 
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the family of constant amplitude A, the function’ ^(x°) is 
obtained as shown in Fig. 3,123. The function lA (sc 0 ) as a single- 
valued function of x°, which is linear in the range 0 £ | x° j £ , 2 . 4 . 

For an input r(t) = 5 sin 0.5t, the computer solution is 
shown in Fig. 3.24.- The amplitude A and the frequency 0 of the 
limit cycle are slowly-varying quantities according to equations 
3 ,54 and the slowly-varying variable x°. , Their average values, 
however, are close to that which can he predicted from the parameter 
plane diagram of Fig. 3.22? i.e., A = 2.8 and O - 4.5 rad/sec. 

■This can be concluded from the diagram (a) of Fig. 3.24. On 
the diagram (b) , the output signal, c (t) is shown Whereby the 
limit cycle is largely attenuated by the block'G-^s) of Fig. 3.21.' 
The low-frequency component in the signal c (t) represents the ■ 
input r (t) - = 5 sin 0.5t at the output of the system. 

Of course, if the input r(t) is not 'present, the system 
will exhibit a limit cycle which can be determined from the inter- 
section of , the M locus x° s= 0 and the £ ** '0 curve on Fig. 3.22 
as x - A sin t, A = 2.6, 0= 4.8. This is checked by the analog 
computer simulation and the obtained solution is shown on Fig. 3.25. 
3 . 6 Conclusion 

The parameter' plane method has been used to indicate ex- 
istence of asymmetrical oscillations in nonlinear control systems. 

A procedure has been developed to determine the oscillations for 
different values of system parameters and input signals. It has 
been shown how a limit cycle can modify the nonlinear character- 
istic for slowly-varying signals. This modification may be of 
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importance when a high-accuracy control system has 'to be de- 
signed in the presence 'of nonlinearities with excessive dead 
zone, hysteresis, backlash, etc. The design technique can be 
directly applied to a large class of plant-adaptive control s'y; 
terns where a sinusoidal signal is used as an identification sig- 
nal. • 

in a future study, the technique may 'be extended to the in- 
vestigation of transient asymmetrical oscillations. Thus,, to 
study how these oscillations are established after 'certain ampli 
tude perturbation, this study should he largely based upon the 
material presented in the following chapter. 

It may also 'be shown Cl6, 17] that the presented analysis. can 
be extended to the case when the signal superimposed on a sinusoid 
is not .only a constant or slowly-varying sinusoid, hut also when . 
the additional signal is described as a Gaussian process, pro- 
vided that the amplitude or standard deviation of the additional 
signal is of no consequence in the analysis. This further gener- 
ates the idea of applying the dual-input describing function. . 

. [15,17 J along with the parameter plane method, ..and investigates 
the case when the input to a nonlinearity of the system' is' a • 
combination of two similar sinusoidal signals. 



3-R1 


REFERENCES 

1. ' MacColl, L.A. .. Fundamental Theory of Servomechanisms , D. Van 

Nostrand 'Co. , Inc., Princeton, N.’J. , 1945; pp'. 78-87. 

2. Lozier, J.C-. , Carrier-Control-led Relay Servos, Elec. . Engr. , 
vol'. '69, Dec.. 1950; pp 1052-1056. 

3. Tsien, H.S;., Engineering Cybernetics . McGraw-Hill Bpok, Co., 
Inc., N. ’Yl , -1954;. pp. 73-82. 

4; Popov, E.P., Self-Sustained Oscillations in the Presence of 
Slowly-varying Input Signals (Russian) , Trans, of the Acad . 
of Science USSR, T. 98, No. 4, 1954. 

5; Pospelov, G.S., Vibrational Linearization of Relay Systems 

(Russian) , Proc. of the' 2nd Nat‘1 Conf. on Automatic Control, 
Ed. by Acad, of Sci, USSR, t. l r , 1955. . 

6., Riabov, B.A., Regime of Self-Sustained Oscillations in Systems- 

- with Asymmetrical Characteristics (Russian), Automatika , Ed. 
by Acad." of Science USSR, ,No. 2. 1956. 

7. Starikova, M.V., Asymmetrical Self-Sustained Oscillations in 
the Presence of External Inputs (Russian) , Automatic Control 
and Computer Techniques , Ed. No. 2; Maschigiz, 1959. 

8. Popov, E.P., and Palitov, I.P., Approximate Methods for Analyz- 
ing Nonlinear Control Systems (book in Russian) , State Press 
for Physics and Mathematical Literature, Moscow, 1960. Avail- 
able in English trans.. Report No. FTD-TT-62-910 from Foreign 
Tech. Div., - AFSC, 'Wright-Patterson Air Force Base, Ohio. 

9. Loeb, J.M., A General Linearizing Process for Nonlinear Contro] 
Systems, in the book. Automatic arid Manual Control , Academic 
Press., Inc., New York, N.Y. , A Tustin, Ed., 1952, pp. 275-283. 

10. Oldenburger, R. , Signal Stabilization of a Control System, 

■ Transactions, American Soc. of Mech~. Enqrs ., New York, N.Y. , 
vol. 79,. No. 8, Nov. 1957,, pp. 1869-1872. 

11. Oldenburger, -R. , and Liu, C.C., Signal Stabilization of a 
Control System, AIEE Trans. , pt. II (Applications and Industry) 
vol. 78, May 1959, pp. 96-100. 

12. Oldenburger, R., and Nakada, T., Signal Stabilization , of Self- 
Oscillatmg Systems, ’ Transactions. Prof, Group on ■Automatic 
control I PE . New York, N. Y., vol. AC76, No; 3, Sept. 1961; 
pp. 319-325. 

13. Sridhar,' R., Signal Stabilization of a Control System with 
Random Inputs, Ph.D. .Thesis, Purdue University, Lafayette, 

Ind. , Jan. ,1960. 


3-R2 


14. 

15. 

16. 
i7. 

"'18. 

19. 

20 . 

21 . 

22 . 

23. 

24. 

25. 

26. 
27. 


Oldenburger, R. , and Shridar, R. , Signal, Stabilization of a 
Control System with Random inputs, AIEE Trans., pt. 11 
(Applications and Industry), vol. No. Nov., 1961, pp260 267. 


West, j.c. , Analytical Technique for Nonlinear Control Syste ms, 
English Univ. Press, London, 1960. 


Boyer,- R.C., Sinusoidal Signal Stabilization, M.S., Thesis, 
Purdue University , Lafayette, Ind., Jan. 1960, 

Gibson, J.E., Nonlinear Automatic Control , McGraw-Hill Book 
Co., Inc., New York, N.Y., 1963. 


Gelb, A., and Van der Velde, W.E., On Limit Cycling Control 
Systems, IRE Trans. On Automatic Control , vol. AC-8, April „ 
1963; pp. 142-157. ‘ , ' 


Atherton, D.P., and Turnbull, G.F., Response of Nonlinear 
Characteristics to Several Inputs and the Use of the Modified 
Nonlinearity Concept in Control Systems, Proc. IEEE, ~ vol. HI 
No. 1, Jan. 1964’; pp. 157-164. 


Turnbull, G.F., Atherton, D.P., and Townsend, J.M., Method for 
the Theoretical Analysis of Relay Feedback Systems, Proc. IEEE , 
vol. 112, No. ‘5, 'May 1965. 


Tsypkin, Y.Z., Theory of Relay cont rol Systems, (book in 
Russian) Gostehizdat, Moscow, 1955. 

McGrath, R.J., and Rideout, V.C., A Simulator Study of -a Two- 
Parameter Adaptive System, IRE Trans, on Automatic Cont rol, 
vol, AC-6, Feb. 1961; pp. 35-42. 


Rajarmau, V.-, Theory of a Two-Parameter Adaptive Control System, 
IRE Trans « on Automatic Control , vol. AC -7, July 1962; pp. 

20 - 26 ^ 


Smith, Kelvin, Adaptive Control Through Sinusoidal Response, 
IRE Trans, on Automatic Control , vol. AC-7, Mar. 1962; pp. 
129-139. 


Smvth, R.K., and Nahi, M.E., Phase and Amplitude Sinusoidal 
Dither Adaptive Control System, Proc. of 196 3. - J t :. Automat ic 
Control Conf .. AICh.E., New York 1963; pp. 302-312. 


Horowitz, I. M. , Comparison of Linear Feedback Systems with 
Self -Oscillating Adaptive Systems, IEEE Trans . on Automatic 
Control, vol. AC-9, No. 4, Oct. 1964; pp. 386-393.- 

Siljak, D.D. , Analysis of Asymmetrical Nonlinear Osillations 
m the Parameter Plane,, IEEE Transactions o n Automatic Control, 
vol. AC-11, NO. 2, April, 1966, pp. 82-87. 



